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FUTURE PUBLICATION OF 
NAVAL RESEARCH LOGISTICS QUARTERLY 


The Office of Naval Research will discontinue publication of the Naval Research Logistics 
Quarterly after distribution of the December 1982 issue. 

The journal will continue to be published, however, and will appear as a nongovernmen¬ 
tal periodical to be published with the cooperation of the Office of Naval Research. 

The new publisher will be John Wiley and Sons, Inc. The Naval Research Logistics Quar¬ 
terly will become a copyrighted journal in the Wiley-lnterscience series, and will be sold on a 
subscription basis by the publisher. The cost for a one year subscription will be $60.00 to indi¬ 
viduals and institutions. Wiley-lnterscience publication will be initiated with the March 1983 
issue, Vol. 30, No. I. 

A substantial professional discount will be offered to members of the Operations Research 
Society of America and The Institute for Management Sciences. Members of ORSA or TIMS 
may subscribe to the journal for $20.00 per year. 

If you are currently receiving the NRLQ by paid subscription through the Superintendent 
of Documents, U.S. Government Printing Office, and your subscription is due to expire after 
the December 1982 issue, a refund for the unfulfilled portion of your subscription will be 
issued to you by the Government Printing Office. The new publisher has agreed to honor 
requests at the GPO rate for unfulfilled portions which were subscribed to previously. The 
Superintendent of Documents will notify subscribers if they are due refunds because of such 
termination, and refund checks will be mailed subsequently. Subscribers may then elect to 
extend their subscriptions at the GPO rate by contacting the publisher and enclosing their 
endorsed refund checks in payment for the issues due subsequent to termination. An order 
form to facilitate this extension and/or renewal by the new publisher is enclosed. 

Distribution through officially authorized Navy lists and SupDocs dissemination to govern¬ 
ment repository library centers will be discontinued after the December 1982 issue. 

The Office of Naval Research regrets that it cannot provide support incident to such distri¬ 
bution and that those activities requiring copies for official purposes will be expected to bear the 
necessary costs directly. 

Due to restrictions governing the use of GPO subscription lists, the publisher will be 
unable to follow-up with solicitations mailed with listings provided by the Superintendent of 
Documents; hence, future promotional efforts may only be through professional and trade 
sources. Because of this some current recipients may not receive anticipated additional notices 
from the publisher, and may not maintain the continuity of their collections if they rely upon 
the receipt of subscription solicitations from the publisher. Thus, all present recipients who are 
interested in receiving the journal are encouraged to initiate communications,with the publisher 
promptly so that their names will not be purged by default. 






It is expected that current editorial policies will be continued in the commercially pub¬ 
lished version. Professor Herbert Solomon of Stanford University will be the new editor in 
charge and will be responsible for the review of manuscripts in the future. New manuscripts 
submitted for consideration for publication may be sent to him immediately (4 copies) at the 
following address: 


Professor Herbert Solomon, Editor 
Naval Research Logistics Quarterly 
Department of Statistics 
Sequoia Hall 
Stanford University 
Stanford, California 94305 

The Office of Naval Research is pleased with the recognition of the contributions to 
research promulgated through the Naval Research Logistics Quarterly. As ONR’s own scientific 
journal it has reflected Navy interests in diverse areas of mathematics research which we 
believed would enrich future work in logistics and systems analysis, and has well served its 
intended purpose of providing a forum for mathematicians and logistics engineers to interact at 
the highest scientific level. The willingness of the publisher and scientific community to per¬ 
petuate the journal on a self-sustaining independent basis is profoundly reassuring that its value 
has been clearly established. With the forthcoming transition to the private sector, the journal 
will be able to better serve the entire scientific community in the area of logistics research. 
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ON SOME STOCHASTIC INEQUALITIES 
INVOLVING MINIMUM OF RANDOM 
VARIABLES* 

Peter Kubat 

The Graduate School of Management 
The University of Rochester 
Rochester, New York 

ABSTRACT 

Let Xi be independent IFR random variables and let T, be independent ex¬ 
ponential random variables such that E[X,] - £lY ( l for all i - 1.2. n. 

Then it is well known that £( min (X,)| > £1 min (K,)). Nevertheless, for 
!<<«• !«/<» 

exponentially distributed Y ,'s and for a decreasing convex function 4>(*> it is 
shown that 


£[<!>( min Uf,))l < £!<$( min (K,»l. 
!</<» I<i<» 


1. INTRODUCTION 

The bounds for the mean lives of series and parallel systems under various assumptions 
on component life distribution were extensively treated in the well-known book by Barlow and 
Proschan [3]. 

In particular, if we denote by X,{Y,) the life length of the Ah component in the first 
(second) series system, then the following proposition holds: 

THEOREM 1: (Barlow and Proschan [3, p. 122]). Let X t (Yj) have continuous distribu¬ 
tion F,(Gj) with the mean n,. Let X x . X„ (Ki . Y n ) be independent and F, < G h 

i - 1, ... , n, i.e., F, is star-shaped with respect to G,[F, < G, iff (1 lx)G~' F,(x) is increasing 

for x ^ 0). Then the mean life of a series system using components with lives JTj. X n is 

greater than the corresponding system mean life using components with lives K,. Y„, that 

is, 

(1.1) £[minU,, ...,*„)]> £lmin( Y . . Y„)l. 

In this note, it will be shown that for £, having Increasing Failure Rate (IFR) and G, 
being exponential we get 

(1.2) £[<fr(mina,. XJ] < £[4>(min(F|. YJ1, 

•This work was supported in part by the Office of Naval Research under grant No. CNA SUB N00014-76-C-OOOI. 
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provided that <!>(•) is decreasing convex. A special case of the inequality (1.2) was already dis¬ 
cussed in [2] where bounds for the discounted mean of the series system life was considered. 

2. BOUNDS FOR THE MEAN OF A FUNCTION OF LIFE OF THE SERIES SYSTEM 

Before we prove the main result the following two lemmata are needed. 

LEMMA 1: Let X] — F h i — 1,2, be two positive continuous r.v.’s such that F\ < F 2 
and £[3^] - E[X 2 \. Then, 

(2.1) £[<&(*,)] ^ £14>(Y 2 )1 
for any convex and bounded function <!>(•). 

PROOF: For i - 1,2 we have 

£[<D(Y,)1- f" <t>(x)dFi(x) 

Jo 1 

--4»(x)F,w[" + f“i/,(x)F,(x)dx 

where we denoted 0 (j r) - 4>'(x) and F}(x) - 1 - F,(x). Since <!>(■) is bounded then 

-4>(jc )F/(x) I — c is finite. Moreover, since ♦(•) is convex then 0 (jc) is increasing and thus 
all the conditions of Lemma 6.4 [3, p. 112] are satisfied and 

f 0 ili(x)F,(x) dx < J" o 0 (jt)£j(x) dx. 

The rest of the proof follows easily. □ 

LEMMA 2: Let Z ~ F z (x) - 1 - e~ Kx and X ~ £,(*) - 1 - e~* x with ElZ] > E[X\. 

Then there exists a unique point, say jt 0 , where the density of X crosses the density of Z and 

the crossing is from above. 

PROOF: £[Z) > £[Y] implies n > \. Solving pe Xjr ° yields x 0 “ 

(p — A) -1 ln(/x/\) >0. □ 

Now we are ready to prove the main theorem. 

THEOREM 2: Let X, ~ F, IFR and Y, — G, with G ( being exponential and £(Y,1 - 
E[Y,\ for i - 1, ... , n. Let <t>(■) be decreasing convex bounded function on [0,<»). Then 

(2.2) £[4>(min *,)] < £[4>(min Y,)]. 

I i 

PROOF: Denote by X - min(^,). Let X — F, then F is IFR and there exists an 

exponential r.v. Z with cumulative distribution function (c.d.f.) G such that £[Jf] - £(Z]. 
Furthermore, F < G [2, p. 107], 

By the Lemma 1 we have £[4>(JT)] < £(4>(Z)], and by Theorem 1 we have 

£tZ] > £lK], where F-min(K 1 . K„). Note that Y is also exponential. If 

£lZ] - E[Y) we are done. 
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If E\Z ] > ElK), then there exists a unique point x 0 where the density of Y crosses the 
density of Z and the crossing is from above. Since 

£(<i>(Z)]« f°° 4»(^)X4? _AJr 
•'0 

and 

£l4>m] - f°° dx 

* , 0 

we get 

(2.3) £[<f»(Z)l -£[«!>< K)] - <t>(x) (\e- k * - ne-**) dx 

which can be written as 

(2.4) J [<J>(x) — d>(xo)l (Xe -At — lie' 11 *) dx + d>(x 0 ) J* Q (Xe -Ajr — /ue _MJr ) dx. 

The second term of (2.4) is, of course, equal to zero and the first term written as, 

J" o 0 — ld>(x) — <J>(x 0 )] (Xe -Ax — fie-**) dx 

+ f l<t>(x) ~ 4>(x fl )l (\e~ Xx - fie~ MX ) dx < 0, 

J *o 

is nonpositive, since for x € [O.xol, <t>(x) - <J>(x 0 ) ^ 0 (<!>(*) is decreasing) and 

Xe -At -< 0 (by Lemma 2). Similarly, for x€ (x 0 ,«>), <t>(x) - <Mx 0 ) < 0 and 
ke~ Xx — fie~ l ‘ x ^ 0. D 

When the function <!>(•) is exponential we obtain the following corollary: 

COROLLARY 1: Let X, — F ( IFR and Y, — G, exponential, for r- 1. n, 

f [X,] - £[ K,] - I /n r Then for d>(x) - e~ ax we get 

2 > 

£[exp{-a(min X,)\\ < ——. 

' « + Tm, 


Finally, let us remark that Theorem 2 will hold even when F] is Increasing Failure Rate 
Average (IFRA), since this condition together with G, being exponential still implies that 
F < G 13, p. 107], 

3. EXAMPLES 

Two simple examples will illustrate the theorem. 

EXAMPLE 1: Consider an n-unit series system, where the unit lifetime X t ~ F r The F, 
is IFR and E[X,\ - 1/m,. When the system fails it is replaced by a new system. The cost of 
the new system at the time t - 0 is C. The expected discounted cost of replacement at the 

time of system failure is CE[e~ aX 1, where X — min(Y|. X„) denotes the life time of the 

system. 
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It follows from Corollary 1 that the upper bound for the discounted replacement cost is 

2>. 

c — - -. 

« + 


Suppose now that F ( is Gamma (X,.a,), a, > 1 , for all i - 1 . it, i.e., the probaoility 

distribution function (p.d.f.) of X, is 

, \ a ^ 

,, , Xy(XyJr) ' -x x 

f(x) - —p^Tj— e ' • x > 0, k, > 0. 

In this case F t is IFR, E[X,\ *= ajk, and thus/tty - X,/a ( . Moreover, 

C£le‘“*) < C/(a/A + 1), 
where A — £ (Xy/a,). 


EXAMPLE 2: Let X — min(Y|. X„). All X t \ are independent, have IFR distribu¬ 

tion Fj and E[X,\ = l//u,. Assume that 4>(x) « (ax + /3) _l , a > 0, /3 > 0. Clearly, <J>(x) is 
decreasing convex and bounded on (0,<»). it follows from Theorem 2 that in this case 




£[(aJlf+ /3)-'l < r 
** 0 

A I 


ax + /3 


Ar' u dx 




where A — £/t,, rf — A/3/a and £ t (rf) denotes the Exponential Integral. The numerical 

i 

values of £,(</) are tabulated in [1]. 
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COMPOUND AVAILABILITY MEASURES 

Laurence A. Baxter 

Department of Applied Mathematics ,nd Statistics 
State University of Sew York at Stony Brook 
Stony Brook, Sew York 

ABSTRACT 

An availability measure is the probability that a two-slate system modeled 
by an alternating renewal process is available at one or more points or intervals. 
The concept of availability measures is extended to formulae for the joint pred¬ 
iction of availability and numbers of breakdowns for repairs) of the system 
during a fixed interval. 


1. INTRODUCTION 

Consider a two-state system, i.e., a machine subject to stochastic failure and repair, and 
suppose that the breakdown/reactivation cycle of the machine can be modeled by means of an 
alternating renewal process (e.g., Cox l4l, Chapter 7). Under this assumption, we can derive a 
variety of useful formulae for predicting the reliability of the system. In particular, we can 
determine the distributions of the Sit), the numbers of failures and repairs in a fixed interval 
Ill, and we can obtain availability measures, probabilities of the form pi/it) - 1 V r € 7j 
where lit) - 1(0) if the system is operating (failed) at t and where Tis an index set compris¬ 
ing a (finite) series of points or intervals 12], The dependence between the Sit) and the /(/) 
precludes the simultaneous prediction of numbers of failures (or repairs) and availability using 
existing theory. In this paper, we present a series of formulae which generalize the availability 
measures, enabling us jointly to predict numbers of breakdowns (or repairs) and availability. 
We call these formulae compound availability measures. 

Availability measures are of use in assessing the likelihood that a repairable machine 
modeled by an alternating renewal process is available at specified times. If the probability falls 
below a certain threshold, arrangements for securing an adequate back-up can be made. The 
distributions of the Sit), on the other hand, are of use in deciding on the allocation of 
sufficient repair facilities for a fixed time interval. The formulae presented in this paper enable 
us to consider these probabilities simultaneously. 

In Section 2 we examine the simplest case and consider probabilities of the form pl/ir) - 
1, Sit) - j] and p{Hu) - 1 V u € [/,/ + x]. Sit) - j }. In Section 3 we consider probabili¬ 
ties of the form p(lit + x) — 1, Sit) - y), and in Section 4 we consider probabilities of the 
form p[lit) — 1, Sit + x) — Sit) - j). Some numerical examples for the Weibull/gamma 
process are given in Section 5. It is first necessary to introduce some notation and briefly to 
review some relevant results. 


VOL. 29, NO. 3. SEPTEMBER 1982 


403 


NAVAL RESEARCH LOGISTICS QUARTERLY 



404 


L A BAX U R 


Let F and G denote the distribution functions of the failure times and repair times, 
respectively, and let their Stieltjes convolution be denoted 

Kit) - F • Git) « f' Fit - u) dGiu ). 

Jo 

Further, the n-fold recursive Stieltjes convolution of K is written K {n) for n - 1,2,3.and 

we define A <ol (r) — HO) if t > (<) 0. Let /, #and k be the densities corresponding to £, G 
and K , respectively. The renewal counting functions of the numbers of failures and repairs in 
(0,rl are denoted N\it) and N 2 it) iN 2 it) and /V 3 (/)), respectively, assuming that there is a 
repair (failure) at time 0. The corresponding renewal functions and renewal densities are 
H,U) - £{/V,(f)) and h,it) - dH,it)/dtii~ 1,2,3). 

The distributions of the N,it) are as follows: 

I 1 - Fit) n-0 

( U) p[N { it) - n\ - | F t _ F . *<»»(,) „ ^ l 

(1.2) p{N 2 it) - n) - K {n) it) - K'^it) 0 

j 1 - Git) n-0 

(1.3) p{N>it) “ «) “ | G . K { "-"it) - G • K [n 'it) n> 1 

(see Barlow and Hunter [1]). The point availability of the two-state system is defined as 
A k it) - p{l k it) - 1} where k - 0(1) if there is a failure (repair) at time 0. The interval avai¬ 
lability is defined as R k ix,t) - p[l k iu) - 1 Vw € [t,t + *]). It can be shown that 

(1.4) A,(t) - Fit) + F * H 2 it) 

(1.5) A 0 it) - Git) - G • H 2 it) 

(1.6) /?|(x,r) « fit + x) + h 2 iu) Fit + x ~ u) du 

(1.7) R 0 (x,r) - f Q hjiu)Fit + x - u) du 

where 

Fit) - 1 - £(r)[2l. 

2. COMPOUND MEASURES OF lit) AND Nit) 

The expressions for the simplest compound availability measures, i.e., p{lit) — 1, 
Nit) - j), are of particular interest as surprisingly simple formulae can be obtained for the 
covariances of the /*(/) and the N t it). 


(2.1) 

p{/,(r)- 1. N { it)~j ) - 

F • K'J'it) 

j>0 

(2.2) 

p|/,(r) - I, N 2 it) -j) - 

F • K {J) it) 

j>0 

(2.3) 

p[l 0 it)~ 1. iV 2 (r) - j) - 

F • G • K fJ) it) 

j>0 

(2.4) 

p\l 0 it)~ 1, A 3 (r)-y}- 

1- 

\ F * G • K {j ~ u it) 

3-0 

J > 1. 
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Observe that /j(/j(r) - 1, IV|(t) = j] - p{/[(f) - 1, N 2 it) « j) for all j. This is to be 
expected in view of the identity /,(r) — N 2 (t) - )V,(r) + 1. 

The corresponding formulae for P' k ix,t,j) = p\I k (u) - 1 V w € [r,r + xl, N,(r) = j] are 


as follows: 




(2.5) 

P\ (x.t.j) = i 

1 

Fit + x) 

f k <J] (u) Fit + x - u) du 
Jo 

j-o 

j > 1 

(2.6) 

P\ ( x.t.j) - 

Fit + x) 

f k ij) (u) Fit + x - u) du 

Jo 

j- 0 

j > 1 

(2.7) 

Pi (x.t.j) - _ 

j" bj(u)F(t + x — it) du 

j > 0 

(2.8) 

Pi (x.t.j) - 

1°, 

[ J" fl bj-\(u) Fit + x ~ u) du 

7-0 

J > 1 

where k (j) it ) 

= dK ij> (t)/dt and 



I f git - u) k {J '(u) du j > 1 

V'>"U> y-o. 


It should be noted that the above formulae, and many of those of Sections 3 and 4, do 
not need to treat the cr j = 0 separately, but this has been done in the interests of clarity. 


The following covariances are readily derived from (2.1)-(2.4) above. 

(2.9) cov (/,(r).Af|(/)} - A, • H 2 (t) = A^(t)H^it) 

(2.10) cov (/,(/).M 2 (/)) = A | • H 2 (t) - A ] U)H 2 (t) 

(2.11) cov { l 0 it).N 2 it)) - A 0 * H 2 it) - A 0 (t)H 2 (t) 

(2.12) cov {/ 0 (/).Af 3 </)} - A 0 * H 2 it) - Aq(i) lH,(r) - lj. 


These may be obtained directly, in which case the identity 
£ytf''’(r) - H 2 • H 2 (t) + H 2 (t) 

j- 1 

will be found helpful, or by exploiting the binary nature of the indicator variable. Thus, for 
example, expression (2.10) follows from the identity 

“<^ 1-^1 - E{N > U)U ‘ U) -" - £l "’ ( " l, ' ( ' 1 - 01 
and the other covariances may be derived similarly. 
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Observe that the expectation of the product of the indicator variable and JV 2 (/) is the 
Stieltjes convolution of their expectations. 

3. COMPOUND MEASURES OF lit + x) AND N(l) 

We now present formulae for probabilities of the form p{l k (t + x) » 1, A'(r) = j}. 

(3.1) />{/,(/ + x) = 1, yV,(r)=y} 

= Fit + *) + f Q fit + y) A 0 ix - y) dy 7 = 0 

= f /r <yl (w) F(r + jc - «) r/w 
0 

+ J 0 J* Q k <J> iu) fit + y - u) A 0 ix ~ y) dudy 
+ J 0 f Q o y _|( m) git + y - w) A\ix - y) dudy j > 1 


\ fit) 


(u)du j ^ 1 
7 = 0 


(3.2) pifit + *) = 1, N 2 U) = 7 } 

= Fit + x) + J Q /(/ + y) ^o(* - y) dy 
+ JqJq fiu) git + y — u)A [ix — y) dudy 7 = 0 

= f k (j} iu)Fit + x - u) du 

Jo 

,/;/>><«)/(, + y - m) A 0 (x-y) dudy 
+ J g J o aj (u)git + y - u ) A fix - y) dudy j > 1 

(3.3) p{l 0 it + x) = 1, N 2 it) = 7 } 

= J* o £,(w) ^(/ + x - u) dudy 

+ J Q J 0 bjiu)fit + y - u) A 0 ix - y) dudy 

+ J* o J" Q k iJ) iu) git + y - u) A\ix - y) dudy 7 > 0 

(3.4) p{/ 0 (r + x) = 1, N } it) = 7 ) 

= / 0 git + y) A t (x - y) dy 7 = 0 

= J" Q bj_\(u)Fit + x~u)du 
+ J 0 J Q bj-\iu) fit + y - u) A 0 ix - y) dudy 

k (J) iu) git + y - u) A t ix - y) dudy j > 1. 

Note that, as would be expected, expressions (3.1)-(3.4) degenerate to (2.1M2.4), 
respectively, at x = 0 and that application of the key renewal theorem shows that 

lim p[l k it + x) — 1, N/it) = 7 } ” — ~r — p[Niit) = 7 } 

x—00 Ml p2 

where mi and M 2 are the means of Fand <7, respectively. 
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The corresponding formulae for p[l k iu) « 1 V u € (r_+ x,t + x + w], A,(r) = y) are 
obtained from (3.l)-(3.4) above on replacing Fit + x) by Fit + x + w) and /4*(x — y ) by 
R k i H,x - >0 as appropriate. Further extensions of these formulae to more general compound 
availability measures of the form 

p{l k it\) - / k it 2 ) = • • ■ = /*(/„) = 1, A,(r) - j) for t < r, < t 2 < ■ • • < t n 
are similarly obtained. 

4. COMPOUND MEASURES OF lit) AND Nit.t + x) 

In this section we consider a different type of compound availability measure, that of the 
form p[l k (t) - 1, Nfit.t + x) = y! where N-it.t + x) * Nft + x) - A,(r) is the number of 
renewals in [t,t + x]. 


(4.1) 

p{I\(t) = 1, N\it,t + x) = j) 



= «,(x.t) 

= f 0 fit + u)p{N 2 ix - u) - y - 1} du 

7-0 


+ Jo X fit — i + u)p[N 2 ix — u) = j — 1} dsdu 

J > 1 

(4.2) 

p[l\it ) = 1. N 2 it,t + x) = y) 

= Rfx.r) + J 0 fit + u) Gix - u) du 



+ X J* h 2 is) fit — s + u) Gix — u) dsdu 
= J" Q fit + u) p[N 2 ix - u) - j] du 

y-0 


+ X X /(' ~ s + M )/,|A 3 (x - «) “ y) dsdu 

j > 1 

(4.3) 

pUoit) — 1. /V 2 (r./ + x) — /) 



— R 0 ix,t) 

y-o 


— X X — s + u) p{iVj(x - u) — j — 1)} dsdu 

J > i 

(4.4) 

p{l 0 il) - 1, Aj(f,f +x)« y! 



- /? 0 (x,r) + J" 0 f g h } is) fit - s + u) Gix - w) rfsrfw 

7 = 0 


— J* o J h 2 is) fit - s + u) p{N}ix - u) - y'} r/st/w 

J > • 


Observe that, as would be expected, (4.1) and (4.2) degenerate to (1.1) and (1.2), respec¬ 
tively, at / - 0 whereas both (4.3) and (4.4) vanish. Note further that on applying the key 
renewal theorem to the above formulae we see that 

lim p{l\it) - l, /V|(r,r + x) - j] - lim p{l 0 it) — 1 , N 2 it,t + x)- j] 

AVix) 7-0 

A J Q *!i(u) p\N 2 (x - u) - j - \)du j ^ I 
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and that 

lim p[I\(t) — 1, N 2 (t,t + x) - y} — lim p[l 0 it) = 1, V 3 (/,/ + x) — y) 

* A J i)t(u) p{N}(x— u) = j) du y^ 0 

where /I - ^i/(mi +M 2 >. 4>(t) = Fit)//i, and 'J'(f) = f g </i(u) du. Thus, we see that the 
asymptotic probabilities are equal to the corresponding probabilities for the equilibrium alternat¬ 
ing renewal process. 

Generalizations to formulae of the form p{l k iu) - lV«f [t,t + jt], N,it + w,t + 
w + x) - j) are readily obtained by obvious modifications to (4.1)-(4.4). 

5. APPLICATIONS 

One important industrial application of the alternating renewal process is to model the 
sequences of failures and repairs of electricity generating plant, such as boilers and turbo¬ 
generators. The author’s experience in analyzing such data shows that, typically, failure times 
are identically distributed, as are repair times, and that these random variables do not appear to 
be dependent. Knowledge of various availability measures is a means of assessing the likeli¬ 
hood that, for example, peak demand will be met over one or more periods. The distributions 
of the renewal counting functions, on the other hand, give some indication of the number of 
breakdowns likely in a fixed interval. These may be used, for example, in deciding the levels 
of stocks of those spare parts which are most commonly used in repair work. The compound 
availability measures enable us to make the two sets of predictions simultaneously. 

As an example, we examine the alternating renewal process with Weibull failure times 
and gamma repair times, the scale parameter being unity in each case, i.e., 

fit) = at " -1 e~'°, g(t) = t v ~ l e _ '/r(r)). 


Values of 

/»{/,(/) - 1. /V,(t) - j\ - pU\it) - 1, iV 2 (r) -y} - F * K (J Ht) 

were evaluated for 0 < t ^ 10 and j — 1,2,3, ... for a variety of values of a and r>. The 
numerical integration was achieved by means of the Cleroux-McConalogue algorithm [3] for 
recursively-defined Stieltjes convolutions using FORTRAN subroutines developed by 
McConalogue [5], (See McConalogue ( 6 ] for further details.) The values of F * K {J) (t) 
obtained are illustrated in Figures 1-4 for the following combinations of shape parameters: 

Figure a 17 

1 1 1 

2 2.5 2.5 

3 4.5 2.5 

4 5 1.25. 
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A COMPARISON OF SEVERAL ESTIMATES OF AVAILABILITY 
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ABSTRACT 

We compare several compering estimales of Ihe availability of a system 
which alternates between two states, ’up* and 'down.* in accordance with an al¬ 
ternating renewal process. Both interval and point estimators are compared 
under several special but representative situations The comparison reaffirms 
the validity and robustness of the log-logistic jackknifed estimates However, 
when the point estimates are compared from the intrinsic criterion of probabili¬ 
ty of concentration, the uniformly minimum variance estimate obtained for the 
Markov model performs very well 


I. INTRODUCTION 

In a system which alternates between two states, "up* or "down." in accordance with an 
alternating renewal process, the long-term point availability (often simply abbreviated as availa¬ 
bility) measures the probability that the process is up at a given instant of time far enough in 
the future. This index measures the long-run expected fraction of time that the system 
operates satisfactorily. It has been shown in reliability textbooks, (see. for example. Barlow and 
Proschan [ID, that under not too restrictive conditions, this index equals the ratio of the mean 
uptime to the sum of the mean uptime and mean downtime Such two-state models are con¬ 
sidered frequently in the literature in problems related to design and maintenance of computer 
systems, communication systems, power plants, etc. System productivity is directly related to 
equipment availability, and, therefore, much attention is currently being placed in the industry 
on the assessment and optimization of plant availability. Any availability improvement program 
has to begin with the estimation of current availability. In this paper, we examine several con¬ 
tending estimators of availability, and study some of their properties. 

In a recent paper, Gaver and Chu [2] have studied the behavior of the jackknife method 
for estimating availability under several different probability models for the underlying distribu¬ 
tions. They have pointed out that frequently in availability studies, it will be difficult to obtain 
clear and unequivocal probability specifications. Their results obtained by Monte Carlo simula¬ 
tion show that the jackknife method is reasonably robust for providing interval estimates and 
valid under several different special but representative probability models. In this paper, we 
carry out an investigation similar in spirit to Gaver and Chu 12). Our point of departure from 
this paper is that we consider an additional estimate of availability (Mazumdar (3D, and we also 
consider the intrinsic properties of the estimales from the point of view of point estimation. 
Our main result can be staled as follows: from the point of view of providing interval esti¬ 
mates, the jackknife estimate of Gaver and Chu appears to be suitably valid and robust, and it 
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performs as well or better than the other contending estimates. However, from the point of 
view of point estimation, when the intrinsic property of probability of concentration (Rao IS)) 
is considered, the uniformly minimum variance unbiased (UMVU) estimator seems to perform 
as well or better. 

2. THE MODEL AND THE ESTIMATES OF AVAILABILITY 


We assume that the times spent by the process in the "up" and "down" state are mutually 
independent, and we have at hand a sequence of In observations giving the successive up and 
downtimes, U\, D\, (J 2 ,0 2 , ... . U„, D„ of the system- Except for one case, we shall make 
the usual assumption that the uptimes £/, have the exponential distribution. The downtimes D ( 
are identically distributed, and in Section 3, we shall assume that their common distribution 
belongs to the exponential, gamma, lognormal or the Weibull family. The parameter of 
interest, availability, is given by 


( 1 ) 


E[U) 

£U/| + E[D\ 


where E[U] and £)Z?) are, respectively, the expectations of U i and D t . 


We consider the following estimators of availability: 


( 2 ) 


(a) the maximum likelihood estimate: 

3 „ V 

" m * e r, . k 


U + D 

where U - £ UJn and D - £ DJn. 


(b) the jackknifed maximum likelihood estimate: 


(3) 

: _ V 

t, . t ; 


V + E 

where 


(3a) 

V - £ Vjn, £- £ EJn, 

/-t i- 1 

(3b) 

V, - nU - (n - 1)Z7_ ( , 

(3c) 

l a _- 

(3d) 

E t - nD - (n - 1 )D_,, 

and 


(3e) 



(4) 


(c) the log-logistic jackknifed estimate (Gaver and Chu [2]) 




e 2 

1 + e 2 
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where 


(4a) 

Z = - T Z, 
n Pt 

(4b) 

Z, = nZ - (n - 1) Z_ 

(4c) 

Z. ; = In [ X U] - In 


U/ J 

(4d) 

Z = In (Z7) - In (Z>). 


(d) the uniformly minimum variance unbiased estimate (Mazumdar, [3]): This estimator 
was derived for the case when both U, and D t are exponentially distributed with un¬ 
known parameters K and fi, respectively. Denote UlD by S. Then the estimator is 
given by 


(5a) 


,(«) 


- (n - 1) I 


V 2,-1 s'<i - s) n -'-' 


1 - (n — 1) £ 


n-l+j 

n~ 2 | 

j } n + j 


if S < 1 


ifS > 1. 


The above reduces to: 


An) 


n- 1 


n — 1 


n — 1 

1 

c 

2 

_ c2 i . 

3 

n 


fl + 1 

O T 

n + 2 

1 


2 


3 


S } ~ ... + (-1)" 


n- 1 
n~l 

I 2 " -2 ! 

| n~ 1 j 

if S' < I 


s n - 


and 


(5b) 


^umvu^) ( 


n — 1 
1 


S’ 1 + 


n - 1 
2 


n + 1 
2 


S' 2 - ... + (-1)" 


n — 1 
n — 1 


2n — 2 
n - l 


if S > t. 


Some examples of -4 umvu are given below. When n - 2, 
1/2S if S < 1 
1 - 1/2S"' if S > r 


(5c) A um »u^) 

When n - 4, 


(5d) 


„(4) 


3 s 3 oJ . 1 
3 „_i . 3 


A C _ . — c2 _i_ C3 

4 S 10 5 + 20 4 


ifS s? 1 


l -? s ' , + ¥ s_3 -2? r3 ifS>1 
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(3) the jackknifed uniformly variance unbiased estimate. 
<6a) 'i*.un.vu<*> - ~ Z »0(«) 

n I-1 

where 

(6b) W,(n) - /i/4 umvu («) - in - 1) A_ umvuJ (n - 1) 

A um v u ( n) is given by (5b), 


n-2 

n-2 

n — 2 

1 

c ^ c 2 | . 

3 

n — 1 

n 

n + 1 

1 

2 

3 


if S-, < 1, 


n — 2 /»— 2 n-2 

1 , 2 , n-2 , , 

,(« - i) - i-sr,' + si, - ... + (-l )"- 2 -r sr,'"- 2 ’ 

n —II Inf 2n+4| 


if S_ ; > 1 


In the terminology related to jackknife estimates, the quantities £,, Z,, If, are referred 
to as pseudovalues. Let 0 be an estimator of the parameter 0 based on a sample of size n. Let 
0, refer to the set of pseudovalues in this situation, / - 1,2. .... n and let 

(7) 0 - - T 0,. 

" ft 

It has been stated (Miller [4]), that under certain general conditions, the statistic 

, ov yfn Co- 0) 


~^-r T (0, - 0) 1 

»-> ft j 


has an approximate t distribution with in - 1) degrees of freedom. Thus, the approximate 
two-sided 100 (1 - a)% confidence interval for 0 will be given by (L a , U a ) where 

(9) V„ - 0 + /, „ £ <*, - 0) 2 | ^/n' /2 
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L a 




" 1/2 / 

£ (0, - 0) 2 / n' n 


where t „ refers to the 

'-7 

dom. 



100% percentile of Student's t with n - 1 degrees of free- 


In a recent paper, Rao [51 has commented on the wisdom of using minimum mean square 
error as the only criterion for obtaining optimum point estimates. He has suggested that 
because of mathematical convenience and intuitive considerations, the minimum mean squared 
error could be used as a postulate to derive estimators, but he has recommended that their 
acceptability should be judged on more intrinsic criteria such as the probability of concentration 
(PC), where this quantity measures 

(10) PC = Pr [lEstimate — True Value! < a) 
for different values of a. 

For the purpose of comparing the performance characteristics of the estimates (a)-(e), we 
choose the following three criteria: (i) estimate probability of coverage of the true availability 
using the formula (9) corresponding to a nominal value a; ii) the estimated mean and variance 
of the confidence interval length; and iii) the estimated probability of concentration. (The 
mean squared error was used in [3] to compare the maximum likelihood and the umvu esti¬ 
mate, and we do not repeat the use of this criterion here ) These estimates are obtained by 
performing Monte Carlo simulation experiments where 1,000 synthetic system realizations were 
observed through a cycle of n consecutive up and downtimes. 

Following Gaver and Chu [21, we assume the following distribution forms for the up and 
down times for the purpose of the simulation study. Although these distribution forms 
represent special cases, they are representative of many different practical situations. 

(A) U is exponentially distributed, E\U] — A" 1 ; D is exponentially distributed, 
E[D\ - n~'\ \U, I, |D,I are mutually independent sequences of independent random 
variables. 

(B) U is exponential with E[L) - A -1 ; D is gamma dislribuled with E\D\ - Var 

ID) - (\fkn)~ 2 . where k and n represent, respectively, the shape and scale parame¬ 
ter of the distribution. The successive up and dowtimes are mutually independent as 
in (A). 

(C) U is exponential with E(U) - A -1 ; D is Weibull distributed with the scale and shape 
parameters y and 5, respectively, given by the following density: 

(11) f D (d) - ; r > 0. 

The successive up and downtimes are independently distributed. 

(D) V is exponential with E(U) — A" 1 ; D is lognormally distributed with parameters n 
and <r, given by the following density: 

(12) f D (d) - — exp[-\/2(\nd - n) 2 /<r 2 ]\ d > 0. 

v2wcra 
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The successive up and downtimes are independently distributed. 

(E) U is represented by a long-tailed It distribution, see Gaver and Chu (2). For a fixed 
value of the tail-stretch parameter h, (0 < h < 1), represent an uptime as 


A 

where X is exponentially distributed with £[,T] = 1. The sequence is one of 
independent random variables, themselves independent of (Dj. The downtimes are 
here assumed independent and exponential. 

In case (A), the ratio U/D is proportional to a statistic having the £ distribution with 
(2 n, In) degrees of freedom. Thus, an exact confidence interval for the availability parameter 
A can be obtained in this situation with the help of the F-tables. The confidence intervals 
using the F-statistic were also computed. 

3. NUMERICAL COMPARISONS 

In this section we compare the different estimators using the three criteria described 
above for the probability models (A) to (E). 1,000 Monte Carlo realizations were used to 
obtain the estimates of these performance measures. As far as practicable, the same random 
numbers were used for the purpose of the comparison. We examined the case where n — 15. 
Table 1 gives the properties of the confidence intervals obtained from (9). Table 2 gives the 
properties of the point estimates using the PC criterion of equation (10). 

TABLE 1 — Simulation Results * (1000 Runs) Comparing the Two-sided 95% Confidence 
Intervals of Carious Estimators of Availability; n = 15, t =* 2.145 
(True Value of Availability = 0.9901). 


Underlying 

Distributions 

Coverage 

(%> 

Average Length 

Variance of 
Length 

A. (exponential. 

F:95 4 

1.65 x 10~ 2 

3.65 x 10- 5 

exponential) 

JK.LL:94.5 

1.87 x 10- 2 

8.00 x 10~ 5 

A - 0.01, fj - 1 

JK.UMVU:92.3 

1.54 x 10~ 2 

4.59 x I0- 5 


JK,MLE:93.2 

1.64 x 10- 2 

5.02 x 10~ 5 

B (exponential. 

F:98.4 

1.67 x I0- 2 

2.52 x 10- 5 

gamma) 

JK,LL:94.7 

1.40 x I0- 2 

3.21 x I0- 5 

A - 0.01, n - 3.0 

JK.UMVU:93.4 

1.29 x 10- 2 

2.64 x Itr- 

k - 1/3 

JK.MLE94.8 

1.38 x I0- 2 

2.89 x 10- 5 

C. (exponential. 

F:98.0 

1.67 x I0‘ 2 

2.65 x 10- 5 

Weibull) 

JK.LL94.3 

1.39 x 10- 2 

3 24 x 10-^ 

A - 0.01, y - 1.13 

JK.UMVU.92.5 

1 29 x I0~ 2 

2.71 x !0- 5 

8 - 2.0 

JK.MLE94.0 

1.38 x I0~ 2 

2.98 x I0- 5 

D (exponential. 

F:92 7 

1 68 x I0' 2 

5.30 x I0- 5 

lognormal) 

JK.LL94.9 

2 28 x 10- 2 

4.37 x 10—» 

A - 0.01,*! - 0.50, 

JK.UMVU:92.5 

1 66 x I0- 2 

8 86 x I0- 5 

(7 - 1.0 

JK.MLE89 8 

1 76 x I0- 2 

9.53 x 10- 5 

E. (long-tailed h. 

F:88 2 

1.77 x JO" 2 

6.74 x |(T 5 

exponential) 

JK.LL93.4 

2.41 x 10- 2 

i .60 x ir* 

A - 0.01, M - 1. 

JK.UMVU 92 7 

1.91 x 10" 2 

8.48 x I0-* 

h - 0.2 

JK.MLE:94.0 

2.01 x I0’ 2 

9.18 x 10- 5 


(*F: F-statistic; JK.LL "Log-Logistic Jackknife;' JK.UMVU: "Jackknife UMVU;” 
JK.MLE: "Jackknife Maximum Likelihood Estimate) 
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TABLE 2 — Simulation Results (1000 Runs) Comparing the Probabilities of 
Concentration for Various Estimates of A vailability; n — 15. 

(PC — Pr ( \Estimate - Availability < «)] 


| Underlying 

Distribution 

Estimate 

Probability of Concentration (PC) 

in % 

a - 0.001 

a - 0 0025 

a = 0 005 

a ~ 01 

A. (exponential. 

‘jmle 

22 9 

53 8 

83 5 

9T4 

exponential) 

•llk.mle 

20 8 

52.7 

86.1 

98.3 

X - 0.01. ,x -= 1 


22 5 

53.9 

860 

98 3 



21.5 

52.6 

85 8 

98 2 


At. n 

224 

534 

83 4 

97.7 

B (exponential. 

‘jmlf 

28 1 

611 

890 

986 

gamma) 

Atm* 

26.6 

61 6 

91 5 

990 

A - 0 01, n = 3.0 

8Uv» 

28.1 

62.1 

91.7 

99 1 

* - 1/3 


26.9 

61 9 

91.3 

99.0 I 


At. 11 

28.1 

61 8 

898 

98.7 

C (exponential. 

'jm k 

26.2 

59.9 

89.1 

98.3 | 

Weibull) 

Atm* 

25.8 

59 1 

91.8 

98 9 1 

A - 001. y = 1.13 

A,*,. 

26.9 

608 

92 4 

99.0 

8 - 2 0 

'jfk.umvu 

26.1 

59 3 

91.6 

98.9 


At. ■ i 

260 

59 4 

90 5 

984 

D (exponential. 

‘imle 

205 

47.1 

82 5 

95.1 

lognormal) 

Atm* 

19.9 

47.4 

83 3 

966 | 

A - 0 01, M - 0 50 

A,*,. 

19.1 

469 

84.6 

96.4 

IT - 1.0 

A 

20.4 

47.4 

836 

966 


A ik. 11 

20.0 

47.2 

806 

94 6 

E. (long-tailed h. 

A»<r 

165 

42 6 

73.5 

92.8 j 

exponential) 

Atm* 

152 

39 2 

70 8 

94.5 

A - 0.01. n - 1 

Amhvu 

17.0 

43.4 

749 

95.1 ] 

6 - 0 2 


155 

40.0 

71 1 

94 1 | 

1 At., 

16.5 

40 7 

72.1 

93.0 | 


Table 1 shows that the log-logistic jackknife estimate provides approximate nominal cover¬ 
age of the true parameter value under a wide variety of situations. In comparison, the jack¬ 
knifed uniformly variance unbiased estimate or the jackknifed maximum likelihood estimate 
does not fare as well. However, when the probability of concentration of the point estimate is 
considered, the uniformly minimum variance unbiased estimate performs very well. In many 
cases, it performs better than the maximum likelihood estimate or the jackknife log-logistic 
estimate. These findings are not surprising in view of the more pronounced locally linear 
characteristic of the log-logistic estimate. 

4. SUMMARY AND CONCLUSIONS 

Several competing estimates of availability of a system which alternates between two 
states — "up" and "down" have been compared. The data is assumed to consist of n complete 
cycles of successive up and downtimes. Monte Carlo simulations carried on several special but 
representative situations indicate the general validity and robustness of the log-logistic jackknife 
estimate of Gaver and Chu. However, when point estimates are considered, from the point of 
view of probability of concentration, the uniformly minimum variance unbiased estimate for 
the Markov case performs very well in these situations. 
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ABSTRACT 

In this paper we consider a simple three-order-statislic asymptotically un¬ 
biased estimator of the Weibull shape parameter c for the case in which all 
three parameters are unknown. Optimal quantiles that minimize the asymptotic 
variance of this estimator. £ are determined and shown to depend only on the 
true (unknown) shape parameter value c and in a rather insensitive way. 
Monte Carlo studies further verified that, in practice where the true shape 
parameter c is unknown, using always c with the optimal quantities that 
correspond to c — 2.0 produces estimates, c*, remarkably close to the theoreti¬ 
cal optimal. A second stage estimation procedure, namely recalculating chased 
on the optimal quantiles corresponding to £*, was not worth the additional 
effort. Benchmark simulation comparisons were also made with the best per¬ 
centile estimator of Zanakis [201 and with a new estimator of Wyckoff. Bain 
and Engelhard) 1181, one that appears to be the best of proposed closed-form 
estimators but uses all sample observations. 

The proposed estimator, c*. should be of interest to practitioners having 
limited resources and to researchers as a starting point for more accurate itera¬ 
tive estimation procedures. Its form is independent of all three Weibull param¬ 
eters and, for not too large sample sizes, it requires the first, last and only one 
other (early) ordered observation. Practical guidelines are provided for choos¬ 
ing the best anticipated estimator of shape for a three-parameter Weibull distri¬ 
bution under different circumstances. 


1. INTRODUCTION 

The cumulative distribution function of a three-parameter Weibull variate AT is given by: 
(1) Fix) - 1 - exp(-((x - a)/b] c ) x ^ a 

•Research performed by this author was supported by the US Office of Naval Research under Contract No. N000I4- 
76-C-0723. 
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where "a" is the threshold, or location parameter, "b " is the scale parameter and V" deter¬ 
mines the shape of the distribution. See Johnson and Kotz, |7J and Mann, Schafer, Sing- 
purwalla [II] for comprehensive information about this distribution and its applications. 

Maximum likelihood estimates (MLEs) of all three parameters have desirable asymptotic- 
properties. However, since no closed-form expression for the MLEs exists, the use of a 
numerical iterative procedure is necessary to solve, for a given sample, the corresponding non¬ 
linear optimization problem. This is a computationally difficult problem for which many algo¬ 
rithms fail to provide solutions [19], Other deficiencies of the maximum-likelihood procedure 
as it applies to this three-parameter distribution are detailed in [16], 

Practitioners having limited time and access to analysts or computing facilities would thus 
prefer using closed-form (analytic) estimators. However, only rather inefficient estimators cf 
this type for the three-parameter Weibull distribution have been available in the past [II], 
Recently, more than seventeen closed-form estimators for the three Weibull parameters were 
compared by Zanakis, and the mcst efficient ones identified [20], 

Of the three Weibull parameters, the shape parameter, c, is particularly important because 
it affects the shape of reliability and hazard rate curves (see [7], [12]). and is the primary cause 
of computational difficulties and estimation errors in MLE problems [10], [21] 

Here, we examine an improved simple percentile estimator, c, of the shape parameter, 
with quantiles determined optimally and independently of the other two Weibull parameters. 

2. A PERCENTILE ESTIMATOR OF SHAPE c WHICH IS INDEPENDENT 
OF THE LOCATION AND SCALE PARAMETERS. 

Let i, < h < ... < i„ be the ordered observations from a random sample of size n from 

( 1 ). 


For the two parameter Weibull (when the location parameter is known or zero) Dubey [3] 
proposed a two percentile (p, < p k ) estimator for the shape parameter c: 

(2) c = ln[ln(l - /?*)/ln(l - p,)}/ln{/*//,) 

with i, = /(^|+i and t k = t\ np ^\np,\ being the greatest integer less than or equal to np r 

He showed (as was shown earlier by Lieblein [9] in investigating the extreme-value distri¬ 
bution) that its asymptotic variance is minimized when 

(3) p, — 0.16731 and p k — 0.97366. 

Murthy and Swartz [13] derived unbiasing factors and hypothesis tests for c. Recently, Bryson 
[1] argued that the asymptotic efficiency of r is not especially sensitive to small deviations from 
the optimal percentiles. 

In a recent paper involving simple estimators for all three Weibull parameters [20], five 
estimators of c were compared. These included the following proposed simple percentile esti¬ 
mator, 

(4) c = 0.5ln[ln(l - />*)/ln(l - /?,)]/ln[(z* - - /,)] 
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where />,, p r p K are three quantiles satisfying 

(5) 0 < p, < p, < p k ^ 1 

( 6 ) —In(1 — p t ) = {ln(l — Pj) • In(1 — />*.)} 1 /2 

and t,. i k (/, < i, < t k ) are the corresponding three ordered -.ample observations. (See also 
Johnson and Kotz [71.) 

This estimator requires only three sample observations and is independent of the other 
two Weibull parameters. Thus, when all three parameters are unknown, estimator c is appeal¬ 
ing if the three percentiles can be selected in some optimal way. It was shown before [201 and 
verified here that using (3) and (4), i.e., 

(7) c 2 = 1.494/1n (<4 - r,)]/(r, - /,) 1 
produces a rather poor behavior. 

The purpose of this paper is to derive optimal (or nearly optimal) quantiles p, and p k (p, 
then obtained from ( 6 )) for c when all three Weibull parameters are unknown and to compare 
the estimators based on these optimal quantiles with certain othcs which also require no com¬ 
plicated iterative procedures. Such simple estimators can be used to obtain final estimates or as 
initial estimates for iterative procedures. 

3. OPTIMAL QUANTILES FOR c 


With the use of a theorem of Mosteller [12] and a lemma of Rao [14], the following ana¬ 
lytic expression for the asymptotic variance of c was obtained: 

12)2 


( 8 ) <t>(p,.p k ,c ) = 


4 nc 2 


ft* ~ Qj 


ft - ft 


R, 

1 2 

1 1 1 

+ 2 

Qr'IQi - Q.) 

or' (q, - ft) ft/-' 

<?, -ft ft - ft 

ftr'fft - ft). 



i i 1 1 

l 

1 , 1 

2 

ft -1 

ft - ft ft - Qi 

ft" 1 

Q, - ft ft - ft 

+ er'<ft - ft> 


Qr'Qr' 


where 


and 


R s = pj (1 — p s ) s = i,j,k 

Q s = [—In (1 - p s )V u s = i.j.k 


so that from < 6 ) 


Qj= (Q,Q,) U2 - 

Thus, for a given sample size n, this asymptotic variance is a function of the true shape param¬ 
eter cand only two of the three quantiles. 
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A pattern search procedure 16], extended by Zanakis [21J for bounded-variable nonlinear 
optimization problems with or without derivatives, was employed to obtain, for selected values 
of the shape parameter c, the quantiles p, and p k that minimize ip,, p k , c)/n. The results are 
summarized in Table 1 along with the function values at the point given by (3). 


TABLE 1 — Optimal Results for c 


True Shape 
Parameter 

c 

Optimal Values 

Non-Optimal* 

rr<(>(pi,p k ,c) 

Asymptotic 

Inefficiency 

4>(p’,p^c)/<t>(p,,p k ,c) 

Pi 

_ Bk _ 

rr<f>(Pi,p k .c) 

|M]n| 

0.0086 

IMi 



2.10 


0.0048 

0.9816 


3.194 

3.11 

1.5 

0.0028 


3.155 

12.314 


2.0 

0.0033 

0.9920 


34.976 

3.85 

2.5 

0.0051 

0.9932 

23.215 

81.286 

3.50 


0.0072 

0.9939 

51.070 

164.374 

3.22 

3.5 

0.0092 

0.9944 

99.545 

300.142 

3.02 


0.0109 

■T > ■ ■ 

176.936 

507.770 

2.87 

4.5 

0.0124 

pjsjgf M 

292.965 

809.210 

2.76 

5.0 

0.0137 

ipMS ■ 

458.762 

1229.437 

2.68 

7.5 

0.0179 

; 1 

2522.945 

6194.795 

2.46 


0.0202 

ESS I 

8314.425 

19586.645 

2.36 


*p'“ 0.16731 and p k — 0.97366 


These results reveal that the asymptotically optimal quantiles (p,,p*) depend on the true 
value of the shape parameter in a rather insensitive way, especially for p k . 

It should be noted that using instead (p,,p k ) from (3), which is optimal for the two- 
parameter Weibull distribution, essentially triples the asymptotic variance of estimator c if all 
three Weibull parameters are unknown. This is more pronounced when the true c — 1.5 to 2.0, 
as the last column of Table l indicates. Note also, from Dubey [3] and Kimball [8], that the 
asymptotic variance of c (given by our Equation (2)) when the location parameter, a, is known 
to be or can be set equal to zero, is 0.92 c 2 /rt. Thus, incorrectly assuming that "a" is nonzero 
makes the estimation of c much more difficult than necessary. 

4. MONTE CARLO INVESTIGATIONS 

In order to gain further insight into the behavior of the percentile estimator c, two Monte 
Carlo investigations were made. The first investigated the sensitivity of c to small deviations 
from the optimal quantiles and whether a simple iterative scheme could improve the perfor¬ 
mance of c. The second investigation compared c to a very recently proposed efficient percen¬ 
tile estimator, c, based on all order statistics [181 and the best percentile estimator, c\ found in 
our previous study [201. 

Sensitivity of l 

Since the asymptotically optimal quantiles depend upon the true (unknown) shape param¬ 
eter c in sucfV an insensitive way, a Monte Carlo study with 1,000 replications was made using 
c, defined by (4), (5), and (6), with p, and p k corresponding to optimal quantile values for 
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c = 2.0 and 3.0. Over ihe range of investigation, i.e., c = 0.5 (0.5)3.5, somewhat better results 
were obtained with c* corresponding to p, = 0.0033, p k = 0.9920 than with one corresponding 
to p, = 0.0072, p k = 0.9939. W e define 

(9) c* = 3.643/In[(r A . - /,)/</, - r,)] 

where t, = /|ooo 33 »l+i. h = 0o.»92n n |+i. and /, = /|oii87 B |+i Note that for most sample sizes of 
interest, r, = r, and t k = t„. 

In this investigation, we also examined the following two estimators: 

£ 0 : defined by (4), (5), ( 6 ) and based on the optimal combination of p, and p k from 
Table 1, although this would not ordinarily be possible since the true c value is not 
known. For this reason, we also examined the following: 

c**: a single iteration modified estimator 

(10) £** = £(/?*, p k ) 

with p*and p* selected in a manner prescribed by the value of c*as given in Table 2. 


Table 2 — Rule for Selection of p, and p k in 
Calculation ofc** = c (p’.pl) 


c* 

pi 

Pi 


0.0086 

0.9746 

0.5 - 1.0 

0.0048 

0.9816 


0.0028 

0.9887 


0.0033 

9 


0.0033 

m 


0.0033 



0.0033 



The schedule shown in Table 2 for selection of combinations of values of p* and p k was 
determined from results involving c 0 and c*. Because of the negative bias of c* for values of 
c ^ 1.5, a correcting term (ranging from 0.1 to 0.7) was also investigated for calculating £** as 
a function of c*. This tended to decrease the bias of c" at the expense of higher mean square 
errors. Early simulation results also revealed that for n < 100, perfect choice of percentiles 
(those corresponding to the true but unknown c) reduces the bias but not necessarily the MSE, 
particularly when the true c is large. This theoretical estimator, c 0 , has optimal asymptotic pro¬ 
perties and is not necessarily good for small n and large c values. 

Comparison with Two Previous Percentile Estimators c' and c 

As a benchmark comparison, the proposed estimators c*and £** were also compared via a 
similar Monte Carlo study of 5,000 replications with the two best percentile estimators, c' and 
c, suggested in earlier studies (20 and 18]. 

The best of the simple estimators for the shape of a three-parameter Weibull distribution 
examined in ( 20 ) was found to be: 
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(11) c' = 2.989/ln((r* - a)/(t,- a)] 
where 

(12) a - (/,/„ - / 2 2 )/(r, + /„ - 2/ 2 > 

was the best estimate of the Weibull location parameter and 
2.989 = In [In (I - 0 . 97366 /ln(l - 0 . 16731 ) 1 . 


WyckofT, Bain and Engelhardt (18] recently extended their earlier estimator c [4] via the 
following three step procedure: 

First determine the shape parameter initial estimator (Cl2.20 in (201) 

(13) c„ = 2.989/ln((/|„„ 4 | + ,—/,)/(/,* <1+1 - /,)) 

with p h p k given by (3), and the location parameter estimator 


< 14 > a = U, - r/« ,/<0 ]/(l - \/n ‘°1 


where t is the sample mean. Note the similarity of estimators (11) and (13). Then the 
WyckofT, Bain, Engelhardt estimator is given by: 


(15) 


nkj 


— T ln(r, — a) + 


£ ln(/ r — a) 

r-s+I 


where s - [0.84n] and the constant k„ is tabulated in (41. 


This estimator uses all order statistics and consequently it is not surprising that it was 
found in (18) to have smaller bias and MSE than the much simpler estimator c\ which is based 
only on three quantiles. However, c' performed better in [20] than the earlier version of c pro¬ 
posed by Engelhardt and Bain (41. 


Monte Carlo results comparing the proposed new estimators c* and c** with the previous 
estimators c' and c are shown in Table 3. The nonoptimal estimator c 2 given by (7) was also 
computed in the same simulation, but MSE values were so large and unpredictable for large 
values of c, that no comparisons were necessary. 


TABLE 3 — Comparison of Average Values and Mean Squared Errors of c',c,c*, and c** 


Average Value 

True c 

MSE 

True c 



0.50 

1.00 

1 50 

2.00 

2.50 

3.00 

3.50 

0 50 

1 00 

1.50 

2.00 

2.50 

300 

3.50 


t' 

0.52 

0.99 

1.38 

1.70 

1.95 

2.18 

2.36 

0 008 

0031 

0.087 

0.227 

0510 

1.031 

1.551 


c 

0.53 

1 04 

1.52 

1.91 

228 

2.60 

2.87 

0 008 

0.035 

0 097 

0191 

0391 

0.705 

1.133 


r* 

0.54 

101 

1.43 

1.81 

2 13 

2.45 

2.80 

0 011 

0038 

0 120 

0.313 

0.673 

1.528 

2.258 



0.50 

0.99 

1.44 

181 

2 13 

2.45 

280 

0008 

0053 

0120 

0 310 

0668 

1 523 

2.271 


?' 

052 

1.01 

1 42 

1.77 

205 

2.30 

248 

0.005 

0.020 

0050 

0 143 

0343 

0.714 

1.144 


c 

0.53 

1.06 

1 55 

201 

236 

2.75 

309 

0.005 

0.021 

0.052 

0.131 

0.236 

0.392 

0.804 


l • 

0.51 

0 98 

1.41 

1.80 

2.16 

2 52 

283 

0.005 

0.021 

0 067 

0.187 

0417 

0851 

1 644 



0 50 

0 99 

142 

1 80 

2 16 

2.52 

283 

0 004 

0025 

0 063 

0183 

0.414 

0.848 

1.642 


?’ 

0 50 

099 

1 42 

1.79 

209 

2 36 

2.58 

0002 

0009 

0.027 

0.087 

0.232 

0.518 

0.991 


c 

0.52 

1 02 

1 52 

1 99 

237 

2 78 

3.13 

0.002 

0.008 

0022 

0059 

0129 

0 240 

0.395 


?• 

0 50 

097 

1.43 

1 87 

2 27 

2 70 

3 11 

0 002 

0011 

0033 

0093 

0 217 

0434 

0823 



051 

1.00 

l 47 

1 87 

227 

2.70 

3 11 

0 002 

0013 

0.026 

0 088 

0215 

0434 

0.823 



Number of replications - 5000 



Number of replications “ 5000 
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Estimator's c** slightly better performance over c* is not considered worthy of the addi¬ 
tional effort required. Therefore, the latter is instead recommended in all cases—especially 
when the true c is large, but with caution if also n $ 30 

All four estimators behave similarly for small values of c. particularly as n increases. 
Comparing the proposed estimator c* with our earlier c' we see that the former has a smaller 
bias in almost all cases. The MSE of c' is clearly smaller than that of c' when n - 100 and 
c > 2.50. In problems with large sample sizes n, c' tends to have almost no variance around 
its expectation; thus, for large n, the MSE of c' consists mostly of the bias squared, with almost 
no contribution from its variance. 

Finally, the bias of c* is about equal to that of c. but the MSE of c is smaller than that of 
c* when c > 2.00. Our results about cand c' are consistent with those in (18). 

5. SOME LITERATURE EXAMPLES 

The estimators examined earlier in this paper, along with some other percentile estimators 
used in 120] were applied to seven literature test problems, with known shape parameter. Prac¬ 
titioners may find the results summarized in Table 4 particularly useful However, these results 
should not be generalized since they represent a static picture to a few examples; something 
like setting your watch at 9:00 and keeping it there: it will be an excellent estimate of time 
twice a day. 


TABLE 4 — Literature Test Problem Results 



Problem 

1 

2 

3 

4 

5 


7 


Source 

15 ] 

15] 

IJ7J 

1)7) 

12) 

115) 

115] 


Sample Size 

40 

40 

100 

100 

100 

100 

100 


a 

10 

20 



0 

1.975 

0 


5 b 

100 

100 



1 

47.072 

47.072 


H c 

2 

3 

2 

2 

12 

1.328 

1.328 


Zanakis c 2 - C8t 

7 697 

2.027 

1 360 

1 577 

0986 

1.546 

1.344* 

3 

Zanakis C12.21t 

2.950 

3.633 

2.001* 

2.1.39 

1.305 

1 436 

1.600 


Zanakis C12.234 

1 801 

2.226 

1.714 

1.747 

1 080 

1,350 

1.377 

i 

Hassanein C25.231 

1 502 

1 933 

1 755 

1.154 

1.170 

1 321* 

1 041 

•f 

Engelhard! 









& Bain C27.23* 

1.657 

1.962 

1.770 

1 225 

1.175* 

1.262 

1.139 

- 

Wyckoff, 








i 

Bain & 









Englehardt c 

2 185 

2.571* 

1.970 

1.870* 

1 280 

1.380 

1,444 


Zanakis-Mann ?* 

2.026* 

1.939 

1712 

1.575 

1.137 

1 231 

1 384 

(p, - 0.0033, p k - 0.9920) 








MLE for c 

2.48 

2.330 

1.783 

1.857 

1 201 + 

1 330+ 

1,767 


•Best percentile estimate for problem. 

+MLE better than best percentile estimate for problem 
^Notation used in [201. 


It is interesting to note that in five out of the seven problems, one or more simple per¬ 
centile estimator was more accurate than the MLEs (iteratively obtained on a computer). As it 
was shown earlier (19), 121], this occurs mostly in problems with small shape parameter values, 
particularly when the sample size is small —a fact of particular interest to practitioners having 
limited resources. 
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6. CONCLUSIONS 

A simple, three-order-statistic estimator c* is developed for estimation of the Weibull 
shape parameter c, when all three parameters are unknown. It was found to compare very well, 
in terms of mean squared error, with an estimator of the same form based on an optimal selec¬ 
tion of order statistics, dependent only upon the unknown true value of c. Compared to our 
previous best simple estimator c' [201, c* has a smaller bias in almost all cases and a smaller 
MSE in problems with large values of n and c. For small values of c. c* and c' are as efficient 
as the new estimator c of Wyckoff, Bain, and Engelhardt (18| which is determined from all 
ordered observations after a single iteration. 

We feel that all three percentile estimators c*, c' and c are useful to practioners and 
researchers, depending on the particular circumstances as suggested in Figure 1. 



Figure I Practical guidelines for choosing an estimalor for the shape of 
a three-parameter Weibull distribution 


The following comments will further clarify and support our suggestions: 

If a computer code is available, one is naturally tempted to use Maximum Likelihood 
Estimation iterative procedures. However, research has shown that: 
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(i) MLEs sometimes do not exist or represent a local optimum [16); 

(ii) The shape parameter is the least accurate of the three Weibull parameters to esti¬ 
mate, causing computational difficulties for MLEs as c gets larger [19); 

(iii) In nonregular cases (c < 2), the simple percentile estimator c' was found to be more 
accurate than the corresponding MLE, particularly when n is small [211; 

(iv) the choice of a percentile estimator like c\ c* or c as a starting point may affect the 
speed of convergence, but not so much the final shape parameter estimate of a good 
iterative MLE procedure [21); and 

(v) If the sample size is too large, determination of MLE or even c is very time 
consuming—a real concern in repetitive situations. 

So, even if a computer is available, one may use a percentile estimator initially, addition¬ 
ally or instead. 

It should be noted that practitioners and students (usually nonstatisticians) often need a 
good estimator of the Weibull shape parameter and do not have the time or access to a com¬ 
puter facility, but only a pocket calculator. In such cases, estimator c— which uses all order 
statistics—will be computationally prohibitive for any but extremely small sample sizes. Instead 
we recommend using c* if n > 50, c' otherwise. The unusual advantage of c* is that it is 
independent of the other Weibull parameters. 

These quick estimators are also well suited for obtaining point and interval estimates for 
the otherwise unattainable optimum solution to large scale integer, combinatorial or nonconvex 
mathematical programming problems [22). A sample of heuristic minimum solutions fits natur¬ 
ally a three parameter Weibull distribution. 

In this paper we also demonstrated how the use of optimization methods can change an 
erratic estimator (c 2 ) into a good estimator (c*), by specifying the appropriate three order 
statistics that minimize the asymptotic variance of this simple estimator. 
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LINEARLY CONSTRAINED PSEUDO-NEWTON METHOD 
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ABSTRACT 

Must Newton-type methods for linearly constrained optimization be either 
of the modified Newton or quasi-Newton variety ’ The contention of this paper 
is that explicity recomputing part of the projected Hessian may be superior to 
both approaches. A computational comparison with MINOS is presented 


1. INTRODUCTION 

In a previous paper [10], we presented a new algorithm for solving the linearly con¬ 
strained nonlinear programming problem: 

Minimize f(x ) 
subject to Ax > b. 

That algorithm is a modified Newton one, i.e., for a problem involving n decision variables, an 
entire matrix of projected second partial derivatives, which may be as large as n x it, must be 
computed at each iteration. Such an approach, which requires order n 2 function evaluations 
and order n } arithmetic operations (multiplications and divisions) at each iteration, may not be 
practical when dealing with larger problems, say n > 20. This paper proposes a closely related 
method, which requires substantially less computational work pet iteration, while retaining the 
desirable features of our previous algorithm. 

The chief novelty of our approach in this paper is the idea of explicitly updating part of 
the second derivative matrix at each iteration. To illustrate the idea, the user chooses a "pipe 
width," n-, of rows to update each iteration, and the "pipe" moves from upper left to lower right 
from iteration to iteration. For example, if the second derivative matrix is 5 x 5 and n « 2, at 
the first iteration the updated elements (indicated by asterisks) would be 

* * 0 0 0 
*•000 
0 0 0 0 9 , 

0 0 0 0 0 

0 0 0 0 0 
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at the second iteration, 

0 0**0 
0 0 * • 0 
.... 0 
.... 0 
0 0 0 0 0 

and at the third iteration 

• 0 0 0 * 

0 0 0 0 * 

0 0 0 0 * 

0 0 0 0 * 


Broadly speaking, Newton-like methods fall into two categories—modified Newton ones 
and quasi-Newton ones. As indicated above, modified Newton methods update an entire 
second derivative matrix at each iteration, while quasi-Newton algorithms iteratively update an 
approximated matrix by adding a matrix of low rank (typically rank 1 or rank 2) to the last such 
approximation. While modified Newton methods generally perform well, they are impractical 
for larger problems. Quasi Newton methods require less work (usually order n function evalua¬ 
tions and order n 2 arithmetic), but are potentially subject to a number of numerical shortcom¬ 
ings, chiefly due to the second derivative approximation process, related to scaling the ele¬ 
ments, restarts, and inability to deal with regions of nonpositive curvature. Our idea is to take 
a middle ground between those two approaches, by letting the user decide how much of the 
second derivative information to recompute each time. Our mathematical analysis shows that 
such an algorithm is a viable idea, and some computational evidence indicates that it may yield 
solutions more rapidly than either equivalent modified or quasi-Newton methods. 

The present method inherits four major strengths from our modified Newton approach. 

First, the theoretical convergence results carry over to the present algorithm. Thus, 
under relatively mild conditions, the algorithm can be shown to converge to a point satisfying 
first order necessary optimality conditions. Under somewhat stronger conditions, the methods 
find a point also satisfying second order necessary conditions. The rate of convergence is also 
shown to be superlinear, that is (Ilx* +I - 3f|j/IU*- x||) — 0 as k — 11 «\ or of order two 
(IU* +1 - jell < c||x* +l - 3c|| 2 for some c), depending on the assumptions made, where (**} 
is the sequence of points generated by the method and 3e is the optimal solution. 

Second, the method is reasonably easy to use. It does not require the user to provide 
analytical derivatives. However, while we have had some success in using this approach on 
problems whose derivatives are not smooth (May, Shocker and Sudharshan [11]), the conver¬ 
gence proofs do assume such properties. Our computational experience indicates that the per¬ 
formance of the method is sensitive to only a few parameters, but will converge dependably 
even with highly nonoptimal parameter values. 

Third, the algorithm can deal directly with areas of nonpositive curvature. When minim¬ 
izing a function, and using a Newton-like method, one would like the second derivative matrix 
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(or its projection, if appropriate) to be positive definite, since then the usual search direction 
generation approach yields a descent (improvement) direction. When the second derivative 
matrix is not positive definite, though, a modified approach is required if we want to be sure to 
obtain a useful search direction. Now the second derivative matrix is stored in a factorized 
form, so it can be efficiently used in the solution of the set of linear equations yielding the 
search direction. In our modified Newton method, we utilized the modified LDL r factorization 
of Gill and Murray [6], where £ is a unit diagonal lower triangular matrix, and D is a diagonal 
matrix with all positive entries. That factorization, when applied to a nonpositive definite 
matrix A , will actually find the factors of A + £, where E is a diagonal matrix with positive 
entries, instead of breaking down, as the usual LDL T factorization would. In our present algo¬ 
rithm, the Gill and Murray LDL T factorization could not be used, since we could not find a 
way of updating it in order n 2 arithmetic when the second derivative matrix is modified by the 
addition of a rank-1 matrix. Instead, we use the factorization of Bunch and Parlett [2], which 
they denote Q T MDM T Q . where Q is a permutation matrix, M unit diagonal lower triangular, 
and D block diagonal, i.e., it may have entries in the first sub- and super-diagonal, with only I 
x 1 and 2x2 blocks allowed. (We will suppress the permutation matrix Q to simplify the 
notation.) Sorenson [17] has shown how to modify the MDM T factorization when a rank 1 
matrix is added to the matrix under consideration. While modified Newton methods often have 
provisions lor dealing with the case of nonpositive definiteness, using special line search tech¬ 
niques, our method appears to be the first one which extends these methodologies to a quasi- 
Newton environment. 

Finally, the algorithm is a "least-constrained" method, in the sense of Lenard [8], which 
means it may find the optimal constraint set, and thus converge, quickly. Since our approach 
uses an e-active constraint strategy, it is desirable to identify those constraints met with equality 
at the optimal solution as quickly as possible. A least constrained method allows for the drop¬ 
ping or adding of several constraints from the active set at each iteration, as opposed to the 
usual simplex-technique related strategy of dropping or adding only one at a time. 

2. THE COORDINATE SYSTEM 

A basic idea in our approach is the representation of the locally feasible region in the 
neighborhood of a point x by a matrix factorization of the active constraint matrix N. That is, 
the columns of N are the normals to the constraint considered active at x. Since we always 
assume that A r is full rank, for an n x u matrix A, we can find an n x n orthogonal matrix Q 
and a u x u upper triangular nonsingular matrix R such that 

R 

QN - 

0 

where 0 is a matrix of zeros. 

Now the last (n — u) rows of Q are all orthogonal to the columns of A, which means that 
they define a set of vectors spanning the linear manifold defined by the active constraints. 
Moving away from x along any such row keeps all the active constraint active. Also, the gen¬ 
eralized inverse of N, N + , is given by N + - [/? _l 0]@. The rows of N + have the property that 
movement away from x along the first row of A + , say, corresponds to dropping the first active 
constraint while retaining all the others. 
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The rows of N*. plus the last (// - u) rows of Q , constitute our coordinate system for 
performing an iteration of the method. Since N changes with any change in the active set, so 
does the coordinate system. The method is thus adaptive, adjusting itself to the local structure 
of the feasible region. This adaptability, though, implies a need for proper second derivative 
matrix maintenance, as outlined in Section 5. 

Note that by the orthogonality of Q, the two sets of rows are orthogonal, but that the 
rows of N + are not necessarily mutually orthogonal. 

3. THE ALGORITHM 

The algorithm requires positive real numbers a./S.y.e.r. and p, with p < I and 
0 2 < p/2p. a and 0 are used for move point acceptability, and p is used for a second-order 
Armijo type test for Newton point accep'ability. An e-active constraint strategy is incorporated 
by letting any column A t of A satisfying (A,) T x - b t < es be considered active, y is of the 
order of machine precision, and is introduced to avoid numerical difficulties, r is required in 
the theoretical development as a lower bound on the negative curvature curve stepsize but 
could be set arbitrarily small. Also required are an initial feasible point, .v n . an initial stepsize 
s° > 0, an initial //-vector or n of -t-l’s and -l's, a symmetric matrix G", and 7r, the number of 
Hessian rows / columns to be updated each iteration. 

Let ,V + and Q be determined by QR factorization of N. Denote the / th column of [/V + ] 7 
by n, and / th column of Q r by </,, that is, 

[/V + ] T = [»| //!... //„] and Q r •= lq\ t/j ... q n \. 

Initialize the current solution point x = ,v°, the current stepsize s - s°, the current direc¬ 
tion sign indicator a = a n . the current matrix of QR approximate second partial derivatives 
G = G°, the sequence index A = 0, the last updated column / = //. and the local variations 
failure indicator / = 0. 

General Iteration 

Step 1: (Determine active constraints.) Set the square of the gradient norm rj 2 = 0. 
Determine N for the current x and s. If the active set has changed since the last 
execution of this step, update Q , /V + , and G. 

Step 2: (Computer second order multiplier approximations.) For each constraint /, 
i - 1,2, ... , w, compute A/, an approximate Karush-Kuhn-Tucker multiplier, by 
evaluating / at two feasible points, a stepsize and a half a stepsize along n h and 
using a three point forward difference approximation. If, for any 0, A/, > 0 and 
neither of the two points evaluated along n t yields an improvement, dropping con¬ 
straint / would be undesirable, so set rf, - 0. For each / such that A f < 0 or one 
of those two points does yield an improvement, compute an approximate second 
partial derivative A 2 /, set d, - lA/ji, and let rj 2 - 17 2 + (A/-) 2 . 

Let x B denote the point with lowest objective function value amongst all those evaluated 
at the step. 
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Step 3: (First and second derivatives for manifold directions explicitly updated). Let q, 
denote the row of Q corresponding to the y'th row of G, which is to be explicitly 
updated this iteration. Then evaluate / at two points, one a feasible stepsize along 
q, away from x and the other a feasible stepsize along — q, away from x. Using a 
central difference formula, compute g t and G jr approximate first and second par¬ 
tial derivatives along q, (the discrepancy in subscripts is due to the fact that we are 
interested in only the last (n - u) rows of Q ). Set 17 2 = rr + g 2 , adding g , to the 
norm of the (projected) gradient. 

Update x B to represent the better of x B from the last step and the lowest function value 
point evaluated at this step. 

Step 4: (First derivatives for manifold directions not updated.) Let q, denote the row of Q 
corresponding to the j th row of G, which is not to be explicitly updated this itera¬ 
tion. Then evaluate / at one point, a feasible stepsize along cr i q l away from x. 
Using the value of G# from the last iteration, and the new function value, com¬ 
pute a second order approximation to g r Add g 2 to i) 2 . 

Update x B if appropriate. 

Step 5: (Compute mixed second partial derivatives.) Consider each row q, as in Step 3. 
Then for every other q k corresponding to the / th row of £7, / < j, evaluate /at a 
feasible point appropriately chosen along q, + <j*, and compute the mixed second 
partial derivative G Jt (this fills in the "pipe" as in Section 1). 

Update x B if appropriate. 

Incorporate all the new second derivatives into the matrix factors of G. 

Note that (n - u)/n passes through this step, with the same A, will explicitly update all 
of G, using — in - u) 2 - (n - u) evaluations. In the worst case (u = 0), then, y n{n - 1) 
evaluations are used, on the average. 

The factors of G may be updated at this point using any rank-one or rank-two formula 
which does not force positive definiteness on G on an infinite number of iterations. 

Updating row/column / of G, for / ^ 1, is equivalent to adding a rank-two update of the 
form wej + e t w r , where w /+ , - 0, is the /th column of the identity matrix, so that 

if a scheme utilizing this special structure was devised, at most n updating passes using 
Sorenson’s method [171 would have to be performed. We have not yet devised a procedure 
exploiting this structure. Our code uses a strategy that is reasonably efficient when 

n < y (n - w), and does two rank-one modifications for each row/column updated—the 

second of which wipes out the undesired results of the first. Using the operation count in [18], 
between 2ir(n— u) 2 + 0(n) and (ll/3)rr(n — u) 2 + 0(n) operations are required for this 
updating. Full details on updating the Bunch-Parlett factorization are given in [17). 

Step 6: (Check accuracy of derivative approximations.) If G is positive semidefinite and 
the gradient norm tj is sufficiently small, stop; x is an acceptable solution. 
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If G is positive definite but the stepsize is too big relative to the gradient norm, i.e.. 
s > (rj/a), set the failure indicator r = 1 and go to Step 12, where the stepsize will be reduced 
with no movement made. Otherwise, scale the Newton search direction d relative to the S'* 
directions. If G is positive definite, go to Step 7, otherwise to Step 8 . 

Step 7: (Eigenanalysis when Hessian approximation is not positive definite.) Determine Z 
and A such that Z is orthogonal, A is diagonal, and D = Z\Z T . Let A w denote 
the most negative eigenvalue of D. If A w = 0, D is positive semidefinite; go to 
Step 8 . Otherwise, a direction of strictly negative curvature exists. Solve 
M T y = Z Q for it. If the gradient norm is too small, go to Step 9. 

Step 8: (Compute direction of line search.) Solve MiD + ZEZ 0 M r b = —g for 8, where 
£ is a diagonal matrix given by £ ;J = (max{ IA /; i , yin - t/)X,y] - A,,! for 
/' = 1,2. n — u, and X = max |A„• I. E is thus the "smallest" diagonal 

I < i < n - u 

additions matrix necessary to force positive definiteness on G. If D is positive 
definite, go to Step 10 and search only along the usual Newton direction. If D is 
not positive definite, consider either the "fixed-up Hessian" direction 
l/VT'r/ + Qlh or the quadratic curve {.v,l.v, = .v + i[i.W*) T d - Pii.el + 
t 2 [iN + ) T d + where Q n is the submatrix consisting of the last in - u) rows 

of 0. Note that (A '*) 7 d - Q',g is just the negative gradient in QR coordinates. 
The second vector, iS' + ) T d + Qly. is guaranteed to be a descent direction for / if 
A w < 0 and the derivative approximations are sufficiently accurate. Our curve 
search is thus similar to that of McCormick [121 and More and Sorenson [14], 
except that we reverse the t and t 2 terms. Our motivation for doing this is that, 
unlike them, we allow a stepsize t to exceed one. Since the gradient direction is 
most useful only close to the current point, we wanted it to be multiplied by t 
rather than t 2 . 

If the curve is to be searched, go to Step 9; if the "fixed-up" Hessian direction, go to Step 

10 . 


Step 9: (Search along a curve.) Use a one variable search method to find a ^ t such 
that the point jc,. satisfies the two term Armijo-type acceptability condition 
fix,.) - fix) < plr'v 2 + it*) 2 y T g + y(/*> 4 .v r C>], trying t - 1 first. If success¬ 
ful, let x B = x,.. In either case, go to S ep 11. 

Step 10:(Regular line search.) Search for a t > 0 such that fix + t{(\*) T d + Qog\) ~ 
fix) < pltid T bf + 8 r g) + y r 2 8 r (7], trying r- 1 first. If successful, redefine 
x B . 

Stepll.df the move point is good enough, reduce the stepsize.) If fix B) - 
fix) > - a 2 /3 2 s 2 , there has not been sufficient function value decrease; go to Step 
12. If fix B ) — fix) < — a 2 fi 2 s 2 , set the local variations failure indicator r — 0. If 
fix b ) - fix) < P 2 y 2 , there has been, in addition, sufficient decrease relative to 
the gradient norm, so choose a new smaller stepsize s'and go to Step 13. Other¬ 
wise, x — x B \ go to Step 1, having moved to x B , but using the same stepsize on 
the next iteration. 
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Step 12:(Failure; try reversing the local variations directions.) Reverse the cr direction 
indicators. If r = 0, this is the first time this step has been encountered with the 
current jcand s; set r — 1 and go to Step 1. Otherwise, if r = 1, let the new step- 
size s' = s/2, set r = 0, x B — x, and go to Step 13. 

Step 13:(Define a new sequence point.) If x ^ x k , k — k + 1, x k — x , s k — s, x — x B . 
and s — s'. Go to Step 1. 

4. CONVERGENCE 

The asymptotic convergence behavior of our method is similar to that of our modified 
Newton method [10] from which it inherits its coordinate system and searches, and our uncon¬ 
strained method [9] from which it obtains its second derivative matrix approximation approach. 
Proofs for the theorems below are thus obtained by extending results from [9] in the manner of 
[10], so that we omit the complete proofs and indicate only the key points used. 

THEOREM 1: (First order convergence.) If/is continuously differentiable on an open 
convex set containing the bounded set [x\A T x > b , fix) ^ /(jv 0 )}, then the method con¬ 
verges to a point satisfying the Karush-Kuhn-Tucker first order necessary optimality conditions. 

PROOF: Since fix) is bounded from below, we must have that {sM is infinite, unless the 
sequence U*} is finite. If so, we can show that the last point found satisfies the optimality 
conditions. Now assume (s*) is infinite. Then [s*] — 0, and, by the acceptability tests at Steps 
6 , 10, and II, we have that the gradient norm is bounded above by a function of the stepsize, 
so that the projected gradient norm also goes to zero. 

The next two results all require that / be twice continuously differentiable, {x k \ converge 
to a unique point x , at which strictly complementary slackness holds, and that there is some 
finite bound on the norm of the largest element of the true projected Hessian at x, which we 
denote H. 

THEOREM 2: (Second order convergence). If, in addition to the above. V’/H satisfies 
a Lipshitz condition in a neighborhood of x, and r is chosen small enough, v satisfies the 
second order necessary optimality condition that H be positive semidefinite. 

PROOF: This result follows from the explicit updating we do of G If H is not positive 
semidefinite, then U has at least one strictly negative eigenvalue. Since {.vM — x, and, by our 
updating, [G*] — H, eventually our curve search direction will be sufficiently close to the 
eigenvector direction associated with that negative eigenvalue, and, if r is small enough, the 
search at Step 9 will be able to be successful infinitely often. Since r is bounded away from 
zero, this would imply that (/(a/)] — —°o, so that H must be positive semidefinite. The need 
for a "sufficiently small r” is due to the fact that the upper limit of the interval over which t* 
may be chosen to satisfy the acceptability test at Step 9 is a-function of the magnitude of the 
most negative eigenvalue of H ', and r must be chosen smaller than that upper limit for the 
proof to follow. 

THEOREM 3: If, in addition to the above, H is sufficiently positive definite in a feasible 
neighborhood of x, the convergence rate of the algorithm is superlinear. If a Lipshitz condition 
also holds, order 2 convergence holds. 
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PROOF: These two results follow from making our approximated Newton line search 
direction computed at Step 8 be sufficiently close to the true direction which would be com¬ 
puted using exact derivatives. The stepsize / = I will then be feasible and acceptable, and. by 
McCormick and Ritter (13, Theorem 2!, the rates claimed follow. 

5. MODIFYING THE SECOND DERIVATIVE MATRIX 

The algorithm maintains an approximation, G, to the second derivative matrix projected 
onto the coordinates defined by the last (n — u > rows of Q. M and D factors can be main¬ 
tained in the lower triangle of G, and, since G is symmetric, it is stored explicitly in the upper 
triangle of the same array. G can thus be used to aid in adjusting the matrix when the active 
set changes. Two cases must be considered—dropping a constraint from the active set, and 
adding one. 

Dropping a constraint from the active set removes a column from R and a row from A’ + , 
and adds a row and column to G. An identity column is added to G, requiring no arithmetic. 
(While this may tend to imbalance G, since it will intermix a direction for which no informa¬ 
tion is available with those for which much better approximations are known, we compensate 
by making the new column the first candidate for explicit updating.) The amount of work 
involved in updating Q and R will depend on the number of constraints deleted and their rela¬ 
tive positions in the active set. Removing a single constraint would require from 0 (last con- 
strain! dropped) to 3un + — u 2 + 0 (n) multiplications, 3un + y «* + 0(n) additions, and 

{u — 1) square roots (first constraint dropped) to update R , and an additional 
3[(u — 2)(w — 1 )/2l + (w — 1)0(1) multiplications/additions for updating Q , using the 
arrangement in [5]. 

Adding a constraint to the active set adds a column to R and a row to V + , deflates G from 
u x u to (u - I) x (</ - I), and modifies Q n multiplications and n additions are required to 
determine the new column of R, and then an (n — u) vector must be reduced using orthogonal 
matrices. Using straightforward multiplication, the (# — a — 1)2x2 Givens' matrices each 
require one square root, 4 multiplications and one addition. Premultiplication of Q involves 
3/i(« - u) - 3n + 0</i - u) multiplications and additions. Updating of G is simplified by 
using those transformations that will zero out all but the last entry in the vector, as in Gill and 
Murray |7], We would like to determine triangular factors, LL T . for G, since then premultipli¬ 
cation by the successive Givens' transformations would change L to lower Hessenberg form, 
and postmultiplication would restore triangularity. 

Now G is recurred as MDM r . with D consisting of I x I and 2x2 blocks, and difficulties 
arise when a 1 x I block is nonpositive or whenever a 2 x 2 block exists, that is, whenever 
there is a nonpositive eigenvalue Fortunately, a fix up of G to be positive definite, so as to 
have LL T factors, is easily done, since whenever G is not positive definite at this step, the com¬ 
plete eigensystem of D has already been determined during the previous execution of Step 7. 
Note that this modification of G would have to be performed once, before each set of con¬ 
straints is added, and it will tend to modify the negative curvature information we try to main¬ 
tain. Again, a least-constrained approach should tend to minimize such modifications. 

The method used follows ideas in Sorenson [17], and is an attempt to replace G with the 
"nearest" positive definite matrix, i.e., G - MDM r is replaced with G - M(D + ZEZ T )M T . 
The exact effect on the second derivative information in G is not clear, but, since we were not 
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able lo derive any scheme for preventing the indefiniteness of G which would be more efficient 
than an entire refactorization, we opted for the fix-up approach. Note that this fix up would be 
exactly what we would have obtained if we had used the Bunch-Parlett factorization and forced 
the projected Hessian approximation to be positive definite. 

The fix-up of G and determination of the LL T factor requires between 
— (n — u) (n — u — 1 ) multiplications, no additions, and (n — u) square roots in the case of 

all 1 x 1 blocks, to (n — u) 2 + 4 (n — u) multiplications, ( n - u) 2 + (n — u — 6) additions, 

and (n — u) square roots in the case of all 2 x 2 blocks. The subsequent straightforward mul¬ 
tiplication of the lower triangular factor and its postmultiplication to regain triangularity 

involves 4(« — u) 2 — (n — u + 2) multiplications and 2(n — u ) 2 - (n — u) additions. It 

should be noted that the reduction of the (« - u) vector 


V| 


v 2 


V n-u 


to 11v11„ by successive premultiplication by (« - u) 2 x 2 Givens' matrices of the form 

... 


could be replaced with a single multiplication by 

the matrix 





Cl 

*t 

0 

0 

0 


s,c 2 

-c,c 2 

*2 

0 

0 

... PjP 2 = 

S1S2C3 

-c,s 2 c 3 

-c 2 c> 

*3 

0 


S|S 2 S 3 f 4 

—C|S 2 S 3 C4 

-c 2 s 3 c 4 

-c 3 c 4 

S4 


which is a special lower Hessenberg matrix in the sense of [5], Thus, savings of approximately 
25% of the work involved in premultiplying the lower triangular factor of G could be saved by 
using Lemma IV of Section 2 of [5]. 

If constraints have been added to the active set, N + increases. Due to the complexity of 
updating V + , especially since several constraints are often added at the same time, N + is 
recomputed at such times. Recomputation of /V + requires ~u 2 n multiplications and 

y u(u - 1 )n additions. 
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6. COMPUTATIONAL RESULTS 

The essential motivating factor in this research was our feeling that a partial explicit 
update of the second derivative matrix might prove to be superior to the usual quasi-Newton 
strategy. We coded our method to see whether numerical experimentation would be consistent 
with that conjecture. (A copy of that code may be obtained from the author.) 

There are actually two different major hypotheses to test here. First, we wanted to see if 
varying ir, the number of rows/columns in the update "pipe," could improve convergence 
behavior. We tested this using different values of it in our new code, and by comparison with 
our older modified Newton method, thus holding the algorithmic structure constant. Second, 
we needed to see whether our new method is competitive with other current software. For this 
comparison, we used a state-of-the-art quasi-Newton method, Murtagh and Saunder's MINOS 
[15] algorithm. On a secondary level, we were curious as to how the inclusion of a low rank 
update might affect the results, so we also ran the test problems using Broyden's rank-1 update 
and the complementary DFP rank-2 update. Note that the latter forces positive definiteness on 
G. 


As with similar numerical experimentation with optimization codes, our results cannot be 
considered an absolute conclusion for either hypothesis. The test results using our methods are 
affected by the parameter settings, particularly the starting stepsize. While we tried to find good 
values for each instance, other parameter settings might have been optimal. The path our algo¬ 
rithm took to the solution may differ with w, so that the results may be an artifact of the start¬ 
ing point. Finally, we used the default parameters in MINOS for runs through it. A more 
experienced MINOS user might have been able to speed its convergence by setting its various 
parameters differently. Note that our stopping criterion, since it is a second order one, is more 
stringent than that of MINOS. 

These caveats aside, our experiments do tend to support our conjectures. Table 1 sum¬ 
marizes the number of function evaluations used until convergence (i.e., agreement with the 
optimal / value to 10 decimal places) for the best configuration of our new method, our 
modified Newton method, and MINOS for six popular test problems. The results are 


TABLE 1 — Comparative Function Evaluation Counts 


Problem and 

Best Form of the New Method 

Our Modified 
Newton 
Method 

MINOS 

r- 

Source 

- 

Update 

Evaluations 

Gauthier [3] 
(Colville #7) 

T 

Rank 2 

111 

187 

496 

Hydrazine Equil¬ 
ibrium [4] 

7 

None 

213 

223 

561 

Rosenbrock (16] 

2 

None 

109 (9] 

97 

192 

Shell Primal [3] 
(Colville #1) 

1 

None 

44 

50 

60 

Water Quality [10] 

2 

None 

271 

480 


Wood (31 
(Colville #4) 

3 

None 

245(9] 

447 
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particularly striking for the Water Quality problem, which is a nonconvex 11 variable model 
derived from a model of the Willmette River. The new method performed about the same as 
our modified Newton one on Hydrazine Equilibrium, Rosenbrock's banana function, and Shell 
Primal, and did much better on the other three. 

To further see the effect of modifying ir, Figure 1 shows, graphically, the number of 
function evaluations required for convergence on the four linearly constrained problems. Shell 
Primal is a relatively easy cubic problem, and Figure 1A shows that not much benefit results 
from either increasing it or adding a low rank update. The best performance occurred for 
it — 0 and no update; since we start with an identity matrix for G this is nothing more than 
steepest descent. 

Gauthier’s problem is a 16 variable quartic. In Figure IB the convex pattern with respect 
to it that was so marked for unconstrained problems [9] appears for the "no update" and "rank- 
2 update" curves. That is, an intermediate value of n is superior to either the modified Newton 
(it = n) or quasi-Newton (n * 0) strategy. The pattern for the rank-1 update is different, con¬ 
founded in part by local variations moves. Since the method moves to the better of the points 
found by line (or curve) search or points used for derivative approximations, fortuitous values 
of the stepsize and coordinate system sometimes lead to fast convergence. 

Hydrazine equilibrium is a tightly constrained 10 variable problem. The feasible region is 
small, and the method spends much time changing constraint sets (the optimum is strictly inte¬ 
rior to the inequality constraints). Again, particularly lucky searches on local variation moves 
had an impact on the convergence. The pattern of Figure 1C is a bit confusing, although the 
modified Newton-like strategy of larger jt seems to be desirable. 

Finally, the performance of the Water Quality problem is illustrated in Figure ID. Since a 
rank-1 or rank-2 update cannot be performed unless the same constraint set is used for two 
consecutive iterations, the constant changes of constraint sets meant that such updates were 
never performed. Here the "U" shaped patterns appear again, lending support to the idea that 
some explicit updating is better than none, but that a complete update is usually not 
worthwhile. This is also a rather hard problem; note that MINOS required 2,292 function 
evaluations to converge. 

Table 2 shows the relative percentage of time used by each part of our algorithm. Our 
code is an experimental one, written by the author, so that a more skillful implementation 
might redistribute these somewhat. Test runs were performed on the University of Pittsburgh’s 
DEC 1099 in the usual batch environment. In this situation, times tend to vary by ±15%. 
The standardized times reported in the last two columns were computed relative to the 28.48 
seconds required for Colville’s [3] timing routine (average of 3 runs). 

The bulk of the time required by the new method was for the numerical derivative 
approximations. This was reassuring, since a basic idea in using a nonderivative method, such 
as ours, is to trade off increased computer time for the human time necessary to compute the deriva¬ 
tives analytically. 

Finally, the standardized time ratio for our new method compares favorably with that of 
MINOS. The percentage reduction in time for our method, relative to MINOS, ranges from 
37-70% for the constrained problems and is about 90% for the unconstrained ones. Considered 
together with the additional time necessary to derive and program the derivatives for a code 
such as MINOS, we can conservatively conclude that our new method is at least competitive 
with a representative, state-of-the-art, quasi-Newton algorithm. 


VOL. 29, NO 3, SEPTEMBER 1982 


NAVAL RESEARCH LOGISTICS QUARTERLY 



























I INI \Rl V CONSIR MM I) I'SI I DllAI WHIN MMII(H) 


441 


TABLE 2 — Timing Results 
Time (sec.) 


Problem 

New Method 

MINOS 

Setup; Finding 
Initial 

Feasible Point 

Evaluations 

and 

Derivatives 

Constraint Handling 
and 

Matrix Updating 

Searches 

Standardized 

Time 

Ratio (total) 

Standardized 

Time 

Ratio (total) 

Gauthier 

28% 

38% 

17% 

17% 

00557 

00885 

Hydrazine 

16% 

75% 

2% 

6% 

0.0331 

0.0572 

Rosen brock 

- 

— 

— 

— 

0.0015(19! 

0.0330 

Shell Primal 

56% 

28% 

16% 

< 1% 

0.0109 

00372 

Water Quality 

14% 

61% 

14% 

12% 

0.1131 

0.1921 

Wood 

- 

- 

- 

- 

0.0066119) 

00509 


7. CONCLUSIONS 

Our thesis, in this paper, is that explicitly updating a few columns of the approximated 
(projected) matrix of second partial derivatives at each iteration may yield convergence 
behavior superior to that of traditional quasi-Newton algorithms. We presented a method for 
linearly constrained optimization, incorporating that idea into our existing algorithmic frame¬ 
work. The user specifies the number of rows/column to be updated at each iteration. Theoret¬ 
ical convergence proofs show that such a method retains the properties of our modified Newton 
algorithm—convergence to a second-order Kuhn-Tucker point, and superlinear rate. Numerical 
experiments indicate that this approach is competitive with a state-of-the-art quasi-Newton 
code, and that a small number of explicitly updated rows/columns (usually 2 or 3) can 
significantly improve the speed of convergence. 
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ABSTRACT 

In this paper, we develop efficient interactive methods for the solution of 
bicriteria nonlinear programming problems. The methods do not require 
trade-off information from the decision maker, pose less cognitive burden and 
converge to the "best compromise solution" fast Two methods, called the 
paired comparison method and comparative trade-off method, are presented 
with examples. A real application of the interactive method to a bicriteria prob¬ 
lem that arose in the planning of the cardiovascular disease control program in 
the U.S. Air Force is also presented. 


INTRODUCTION 

Bicriteria mathematical programming (BCMP) problems as a first generalization of single 
criterion nonlinear programs have been studied by many researchers [I, 2, 3, 4, 7, 13, 17J. 
The small dimensionality of the problem permits a visual presentation of the "payoff set," which 
partly explains such widespread interest. Further, when we restrict ourselves to "efficient solu¬ 
tions," the complementary roles of the two criteria yields additional computational power. 

The BCMP problem may be stated as follows: 

(1) VMAX {/, (x). f 2 (x)] 

Subject to £,(*) ^ 0 i - 1. m 

where xis an n dimensional vector of decision variables, ,/j and J 2 are real valued criterion func¬ 
tions and ft's represent a set of nonlinear constraints. Let S - (xlft(x) < 0l denote the feasi¬ 
ble region. An optimal solution of BCMP will be taken to be a "best compromise solution" that 
maximizes the preferences of a Decision Maker (DM); that is, a feasible solution that solves 
the following program. 
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(2) Max U[j\(x), f 2 (x)] 

Subject to x € S 

where U is the DM’s preference function or value function defined over the criterion values 
< / 1 . J 2 ) such that for any two alternatives X| and x 2 , t/l/i (a-j ), y‘ 2 (a - ( )] > [/ ( (jc 2 ), /j(xi)l. if 

the DM prefers X| to x 2 . 

In one of the earlier papers addressed to the bicriteria problem. Geoffrion (7] considered 
the "scalar maximization approach," i.e., maximizing a convex combination of the two criteria. 
He showed how such problems can be numerically solved using parametric programming algo¬ 
rithms. Bacopoulas and Singer [21 on the other hand, used the "constrained criteria approach.” 
i.e., maximizing one criterion keeping the other criterion as a constraint. They showed how all 
the efficient solutions can be generated by parametrically varying the level of the constrained 
criterion over a particular interval. Some further refinements and alternate proofs appear in 
Gearhart [61 and Benson [31. Pasternak and Passy [13], addressed the special case when the 
variables are zero-one and devised a special enumerative algorithm for generating the efficient 
solutions. Choo and Atkins [4] address the case when the criteria are linear fractional func¬ 
tions. Adulbhan [I] considers the special case when the criteria are linear. In all these 
methods, the emphasis is on generating all the efficient solutions of BCMP. But Walker [17] 
considers the "interactive” approach and uses the Generalized Lagrange Multiplier method of 
Everett [5], Walker’s method is slow in convergence and seeks difficult interaction from the 
DM. 


In this paper we will consider interactive approaches to the solution of BCMP where the 
preference function U is only implicitly known. Such an approach in the more general context 
of multiple criteria problems have been addressed by Geoffrion, Dyer and Feinberg [8] and 
Zionis [191. However, Wallenius [18] has observed that such methods pose considerable cogni¬ 
tive burden on the part of the DM; also the procedures converge to the optimal solution rather 
slowly. In this paper, we develop interactive procedures that pose less cognitive burden to the 
DM and converge to the optimal solution fast. 

Section 2 presents the principal results that form the basis of the interactive methods 
developed for BCMP. Section 3 discusses one of the interactive approaches called the Paired 
Comparison Method. Section 4 discusses another approach called the Comparative Trade off 
Method. Section 5 discusses a real application of the interactive method to a bicriteria problem 
that arose in the control of cardiovascular disease in the U S. Air Force. 

2. PRINCIPAL RESULTS 

In this section, we will establish a number of results that form the basis of the new 
interactive procedures which are develolped in the later sections. We will assume throughout 
that the feasible region S is a compact convex set; the criterion functions f\ and f 2 are 
differentiable concave functions and the preference function U is differentiable, increasing and 
strongly quasiconcave defined in the following manner. 

DEFINITION 1: /: S —=> R\ where S is a nonempty convex set in R" is said to be 
strongly quasiconcave on S, if for each X|, x 2 € 5 with x, ^ x 2 and for every X € (0, 1), we 
have 

(3) /[Xx| + (1 - \)x 2 ] > Min [/(x,), /(x 2 )]. 
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The convexity assumptions on criteria and feasible region are made only to allow the 
methods of convex programming to be applicable. The quasiconcavity assumption of L' is well 
supported by preference theory (see Intriligator [10)). 

Use of the concavity assumptions guarantee the attainment of maxima and quasiconcavity 
ensures us that local optima are also global optima. In addition it is sufficient to restrict 
ourselves to efficient solutions defined below as candidates for the optima. 

DEFINITION 2: A solution x° € 5 is said to be efficient to BCMP. if /, (.v) > /, <.x‘') for 
some x € S => / 2 (x°> > / 2 (x). 

We will now develop procedures for generating the efficient solutions and then a method 
of "efficiently" searching among the efficient solutions for the "optimal" solution. 

(4) Define the payoff set Y - |y|/(x) = y for x 6 S). 

(5) Let A = {>’|/(x) ^ v for x € S). Obviously, Y C A. 

LEMMA 1: If 5 is a convex set and /is concave, then A will be a convex set [12], 

Consider the following single objective nonlinear programs: 

PI: Max/,(x) P2: Max / 2 (x) 

Subject to: x € S Subject to: x € S 

Let the optimal values of (PI) and (P2) be v*and w* respectively. Now consider the fol¬ 
lowing mathematical programs: 

Py. Max/ 2 (x) Q w Max,/j(x) 

Subject to: x € S Subject to: x € S 

/|(x) ^ v / 2 (x) > H' 

Let x* solve Q w for w = w*. Then ,/j(x 2 ) is the minimum achievable value for / 
without sacrificing any achievement on / 2 , while v*is the maximum value of ,/j at the expense 
of / 2 . Hence, the range of achievable values for ./j, denoted by v, is given by (,/j(x 2 ), v*] and 
the best compromise value for f\ lies in this range. 

THEOREM 1: In the optimal solution to the program where v € (/j(x 2 ), v*l, the 
constraint ,/j (x) / v will be a binding constraint, i.e., if x solves P- then /|(x) » v. 

PROOF: Assume the contrary, i.e , 

(6) /,(x) > v ^ /,(xj) 

w* being the absolute maximum of f 2 over 5, 

(7) w* - / 2 (xJ) > / 2 (x) 
we will consider the two cases 

(8) (i) Let / 2 (xJ) - / 2 (x) 

(6) and (7) taken together contradict the fact that xj solves Qt. 
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(9) (ii) Let / 2 (x 2 ) > / 2 (.v) 

consider the line segment lx x*l which lies in the convex set S. Over this line segment, the 
concave functions j\ and f 2 are unimodal. Because of “strict" inequality in (6) and (9). there 
exists a feasible solution x' on the line segment (x, xj] such that 

(10) /,(*> > /,<*') > v 
/j(x) < / 2 (x’> < / 2 <xJ). 

In equality (10) indicates that x’ is a feasible solution to £ s with / 2 (.v') > / 2 (x) contradicting 
the optimality of x. 

Hence, f\(x) = v. 

THEOREM 2: x € S is efficient if and only if x solves P- where v € [j\ (,vf), v*l. 

PROOF: 

Sufficiency: Let xsolve £- where v € [/|(x*),v*). We will show that .vis efficient. Con¬ 
sider an ,v € S such that 

/i(.v) > /,(.v) = v. 

It is clear that x is feasible to P- But x cannot be optimal to P x : as per theorem I every 
optimal solution x° satifies /|(x°) - v. Hence, / 2 (x) < ,/ 2 (x). Consider an x 6 5such that 

/ 2 (x) > / 2 (x). 

Since x is optimal to P -, xcannot be feasible to P-. Hence, / ( (x) < J\(x). Consequently, v 
must be efficient. 

Necessity: Let £ be the set of all efficient solutions. Suppose x € £. Obviously f t (x) ^ 
/,(xp as, otherwise, x* will dominate x. Let f\(x) - v. Assume that xdoes not solve P- but 
some other x' € 5 solves £-. By theorem 1, /|(x') — v » /|(x) since x' is optimal for P-. 
/ 2 (x') ^ / 2 (x): but x € £ —> / 2 (x') < / 2 <x) and hence / 2 (x') - / 2 (x). Consequently, x 
solves 

Theorem 2 enables us to generate the entire set of efficient solutions by parametrically 
solving P-. By theorem 1, the generated solutions will have specific levels of attainment of /,. 
namely v. The next theorem would give some unimodality property enabling us to make an 
efficient parametrization. 

Consider the problem where we maximize U over the pay-off set Y, keeping J\ fixed at v; 
i.e.. 


P3: Max U(\,f 2 ) 

<v./ 2 )0- 

where (v,/ 2 ) € Kis a subset of Tdefined by {(v,/ 2 )|/ 2 - / 2 (x) for x € Sand /|(x) - v| 
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By Theorem 1 and the fact that U is increasing in f 2 . it is easy to see that an optimal solu¬ 
tion to (P3) is also an optimal solution to program (/\) where v 6 l./|(.v*). v*]. We now estab¬ 
lish the unimodality of a'(v), defined below: 

THEOREM 3: Let Y' be a convex subset of R 2 and U(J\J 2 ) is a strongly quasiconcave 
increasing function of J\ and f 2 defined over then 

,?(v) = Max i/(v./',) 

k./jK r 

is strongly quasiconcave, i.e.. unimodal in v. 

PROOF: Let 

#(v,) = Max LHv-,./'•>) = 6/<v,.w 3 > (say) 

u ,./,>«> ■'* 

s(v 4 ) = Max t/(v 4 ./,) = 6 (v 4 ,h- 4 > (say). 
u 4 /,ier 

Assuming boundedness of functions the maxima are attained; by convexity of > and strong 
quasiconcavity of U they are attained uniquely. 

Note that (v 3 ,w 3 ) € Y' and (v 4 .w 4 ) € V by convexity of Y'. 

[Av 3 + (1 - X )v 4 , Xh - 3 + (I - A )m' 4 ] also € Y, for 0 < A ^ 1. 
g [v 3 + (1 - A)v 4 ] = Max £/[Av, -f (I — A)v 4 ./ 3 ) 

Iav j + M-x ) 

(By convexity of Y ) > t/[Av, + (1 - A)v 4 .Aw 3 + (I - A)w 4 l 

(By strong quasi- > Min (f/(v 3 ,w ( ), U (v 4 ,h 4 )] 
concavity of 60 = Min (jc<v 3 ), #(v 4 )] 

Hence, g(\) is unimodal in v. 

In the proof of the above theorem, we have heavily exploited the convexity ol > When 
we want to apply the result to Y, the pay-off set, we find that the assumption of convexity of 
pay-off set may be overly restrictive. For example our assumptions of convexity ol 5 and con¬ 
cavity of f/s are not sufficient to guarantee convexity of Y. as the following simple example 
would illustrate. 

EXAMPLE: 5 - x > 0: /= (f,./ 2 ) = (x. -log x) Obviously, 5 is a convex set and s 
are concave. But the pay-off set Y is the graph or - log .v which obviously is not a convex set 
112 ). 

However, we will presently establish the fact that we can extend the pay-off set f to the set .4 
defined in Equation (5) without affecting the results. Lemma I established the convexity of 
the set A. Hence identifying A with Y\ the result of theorem 3 holds even when the pay-off 
set is nonconvex. 
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We will establish in the next theorem, the equivalence of solutions to the following two 
problems. 

P4. Max U ia.f 2 ) P5: Max Uia,f 2 ) 

(a.fjtiY ' («./,)€ 4 

THEOREM 4: The optimal solutions to iPA) and iP 5) are identical, i.e., if Y„ be the set 
of optimal solutions to (P4) and A m be the set of optimal solutions to (P5), then Y m — A m . 

PROOF: 

(a) We will first prove that every element of A m is an element of Y m i.e., A m D Y„. 

Let y* = (yf.y?) maximize P5 for a = yf, i.e., y* € A m . We have to show that 
y* 6 Y and hence y* € Y m . Since A m C A , by definition 3 and x*, such that 
fix*) — (./f,./?) ^ y* Note that y* € Y is f(x*) — y*. Assume that 

fix*) > y*, i.e., 

(ID /?> yf. /? > y* (say). 

Since UiJ\.f 2 ) is increasing in f\ and / 2 , (11) gives 

(12) £/<yf,y?) < Uiytf, ?). 

But from the fact that y* maximizes U(yf,f 2 ) over (yf,/ 2 ) 6 A we get 

(13) (/(yf.y?) > UiyUfi). 

(13) contradicts (12); hence our assumption (11) must be invalid and fix*) - y*. By 
definition of the set Y. y* € Y. Being a maximizer of (PS) for a = yf, y* 6 Y m and hence 
A m 3 Y m . 

(b) Assume y* € Y m . We will show that y* € A m as well or Y m C A m . Since 
Y„, C Y c A, we gel y* € A. Assume y* is not optimal to (P5). Hence 3 
another y'= (y*,y 2 ) such that 

Uiy',v 2 ) •» Max (/(yf,/jl 
a 

(14) (/(yf.y?) ^ (/(yf,/ 2 ) for any (yf,/ 2 ) € 4 
and in particular > (2 (yf.y?). 

Now (yf,y 2 ) € A m ; by part (a) of the proof we will have (y* Y 2 ) € Y. By optimality of 
(yf.y?) for (P4) 

(15) U (yf,y 2 ) < U (yf.y*) 

(14) and (15) imply that 

Uiy*.y 2 ) = V (yf.y?) 

or in other words (yf.y?) optimizes (PS) as well. Hence 
y* 6 A m or Y m C A„. 
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THEOREM 5: With the assumptions mentioned earlier the following relationship is true: 

Max t/(v,/'i) - tr(v) = Maxt/tv,/,) 
i 

Subject to: /j ( v) > v 
.v € S 

PROOF: Follows directly from Theorems 3 and 4. 

3. PAIRED COMPARISON METHOD 

Using the results of Section 2 we will develop a procedure needing onl> a /wired compari¬ 
son from the DM, i.e., given any two feasible solutions and their outcomes. sa> >' 1 ‘ and 
the DM only has to specify whether 

(16) >’ (|1 > >’ UI or y <2) > >•'" or both 

where denotes "preferred to." 

Interaction of this form is presumably much less demanding than asking the DM to 
specify his local tradeoffs. We will presently indicate how (16) can be used to eliminate a por¬ 
tion of the efficient set and progressively converge to the "best compromise solution." 

Based on the theorems given in Section 2, the BCMP problem has been reduced to deter¬ 
mining the maximum of g(v) where v belongs to the interval (v ;> v*l where v ; = ,/j(.v*>. How¬ 
ever, jg(v) is not known explicitly since U is not known. But. using a search technique which 
requires only functional comparison and not function values, we can still solve the BCMP prob¬ 
lem, using the following region elimination concept: 

For a unimodal function jg(v) defined over a finite interval (v,.v*). let v 4 and v fl be two 
points in the interval such that v 4 < v fl . Then, g(v 4 ) < g(v s ) implies that the maximum of 
g(\) will not lie in the interval (v,,v 4 ). On the other hand. g(v A ) > g(v s ) implies that the 
maximum will not lie in the interval (v fll v*). 

General Steps of the Algorithm 

Step 1: Solve PI: Max ./j(.v), subject to v € S. 

Set max /j (x) = v* 

Solve P2: Max f 2 (x), subject to .v € S. 

Set max ./' 2 fx) = w* 

Step 2: Solve Ql- Max ./j(x), subject to * € Sand f 2 (x) ^ tv*. 

Set the maximum of ./j v,. Now the optimal value 
of ./j to the BCMP lies between v, and v*. 

Step 3: Choose two values, v A and v B , such that v, < v A < v fl < v*. Solve the problem 
P v : Maximize ./‘ 2 (jc), subject to x € S and J\(x ) ^ v for v = v A and v = v g ; let 
g(v) — Max / 2 (jv) for P v . 
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Step 4: Let y ni = (v^.gfvy), y (2) = [v fl ,g(v B )]. Given V" and y'*’, the DM is asked to 
specify whether y {|) > (preferred to) y <2 \ or y 12> ^ y m , or indifferent. Usine the 
DM’s response, a portion of the efficient set can be eliminated. 

If y" 1 > y (2) , then set \ u — v B ; go to step 5. 

If v <2> >■ y n , then set v, = v 4 ; go to step 5. 

Step 5: If |v, — vj < « (chosen small number) stop; otherwise return to step 3. 

Figure 1 illustrates the paired comparison method. Here, y 121 > y" 1 since U(y <2 ') > 
U (y 111 ). Hence, the maximum of U cannot lie in the interval [v/.v j. Eliminating the region 
(v,,v 4 ), the interval of uncertainty where the optimum lies reduces to (v^.v*). If Golden Sec¬ 
tion Section Search is used to generate the two new points between v 4 and v*one of the new 
points will turn out to be v fi . With Golden Section Search, we will be able to bracket the best 
compromise solution y° to less than 10% of the original interval on v with just 5 paired compar¬ 
isons from the DM; in 10 paired comparisons, the optimal solution will be within 1% of the ori¬ 
ginal interval of uncertainty [v ( ,v*]. Thus, for a specified level of uncertainty we have a finite, 
rapidly convergent procedure using only paired comparison of two efficient solutions. 
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Numerical Example 

Consider the following bicriteria program [71 with 

/i( x) = 32 —40 x 2 + 23 Xj — 7 x 4 

f 2 (x) = 32 — 10 A| — 4 x 2 + 7 X) — 7 jr 4 

and the feasible region 5 is specified by the following linear equality constraints in detached 
coefficient form in addition to the usual nonnegativity restrictions: 


Xl 

x 2 

*3 

*4 

*5 

*6 

_* 

1 

0 

-23/50 

27/50 

0 

0 


0 

1 

-65/100 

35/100 

0 

0 


0 

0 

1/25 

1/25 

1 

0 

16/25 

0 

0 

15/50 

-35/50 

0 

1 



In order to simulate the interaction process, we will assume that the DM’s implicit prefer¬ 
ence function (which guides him in the selection process) is the following quasiconcave func¬ 
tion 

uiA.fi>-/ Pa. 

The true optimal solution is x* = (jC|,jc 2 ,a 3 ,a 4 ,X 5 ,a 6 ) 

- (1.2857, 0, 0.7895, 3.1805, 0.4812, 2.7895) with ft - (./?, ft) 

« (27.895, 2.406) and U(ft, ft) - 22.13. 

Solution Using Paired Comparison Method 

Step 1. Solving (PI) and (P2) we get v* = 60 and w* « 3.2 
Hence v u = 60 

Step 2: Solving ((?*) yields v, - 8 
Iteration 1 

Step 3: Using Golden Section ratios 0.618 and 0.382, v 4 = 27.864 and v s = 40.136. Solve 
the problem P v for v - v 4 and v fl . g(v^) = 2.407 and g(v fl ) = 1.707 

Step 4: Interact with the DM and seek his paired comparison between y (u - (27.864, 
2.407) and y l2) — (40.136, 1.707). Assuming that DM is guided by the implicit 
preference function U(f\.f 2 ) - f\ /} f 2 , £/Ly <n l — 22.12 and U[y (2) l — 20.01. By 
the property of the preference function, / n > y Q) if t/[> ,<l, l > (/[>’ <2) ]. Hence, 
y u would be preferred to y <2> Consequently, v u will be updated to v s — 40.136, 
thereby eliminating a portion of the efficient frontier. 

Step 5: Assuming that the interval of uncertainty is to be reduced to 10%, the search has 
to be continued further on the reduced efficient frontier. For brevity of presenta¬ 
tion the results of further calculations are summarized in Table 1. 
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TABLE 1 — Solution of Bicriteria Problem Using 
Paired Comparison Method 


Iteration # 

v t 

V U 

V 4 

v 8 




lV :, i 

Remaining 
Interval 
("/.i of Orig¬ 
inal) 

0 

8 

60 

- 

- 

- 


- 

- 

100 

1 

8 

40.136 

27.864 

40.136 

(27.864. 

2.407) 

<40.136. 

1.707) 

22 12 

20.01 

61.8 

2 

20.276 

40 136 

20.276 

27 864 

(20.276. 

2.841) 

(27.864, 

2.407) 

21.12 

22.12 

38.2 

3 

20.276 

32.549 

27.864 

32.549 

(27.864. 

2.407) 

(32.549. 

2.140) 

22.12 

21.82 

23.6 

4 

24.964 

32.549 

24.964 

27.864 

(24.964, 

2.573) 

<27 864, 
2.407) 

21.97 

22.12 

14.6 

5 

24.964 

29.652 

27 864 

29652 

(27.864, 

2.407) 

(29,652. 

2.306) 

22.12 

22.09 

9.0 


4. COMPARATIVE TRADEOFF METHOD 

The local trade-off at a feasible point x (l S, between criteria 2 and 1 is defined to be 



evaluated at x. 

T 2 \ measures the loss in criterion f 2 from the current level the DM is prepared to trade¬ 
off for unit gain in criterion f\. 

If the DM can provide local trade-offs then the problem can be solved interactively using 
Geoffrion-Dyer-Feinberg [8] type approaches. However, it has been repeatedly observed in 
practice that precise tradeoffs are difficult to provide. However, the DM can provide "impre¬ 
cise” tradeoffs much more easily, e.g.. "interval estimates" in place of "point estimates." For the 
bicriteria case, an estimate of the local tradeoff that is easier than the interval estimate is the 
"comparative estimate." The DM is provided with a number and will be asked to respond 
whether or not he would be prepared to tradeoff more than, less than, or equal to the specified 
number in criterion f 2 for a unit gain in criterion f v Based on the DM’s response, a portion of 
the efficient frontier can be eliminated. We will presently indicate how the results of Section 2 
can be used to provide the basis for the "comparative trade-off method." 

At any feasible point x € S, let X 2) be the perturbation in f 2 from its current level per 
unit perturbation in f x . Limiting ourselves only to efficient solutions, we find that Theorem 2 
provides us a method of generating efficient solutions; also Theorem 1 specifies that the con¬ 
straint /|(x) > v in program P- will be binding in every optimal solution. Hence, x 2) is 
obtainable as the negative of the Lagrange multiplier associated with the constraint f\(x) > v. 
Note that X 2) is part of the solution output of most mathematical programming algorithms 
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X 2t has the following interpretation: k 2 \ measures the loss in ,/j from its current level that 
must be suffered per unit gain in /j. On the other hand T 2 \ measures the loss in f 2 the DM is 
prepared to trade-off per unit gain in /,. Naturally, these two quantities must be related to the 
optimality of BCMP and the precise relation is the subject of our next theorem. 

THEOREM 6: 


Let x* be the optimal solution to BCMP with fix ') = (/*, /*). Let x be any efficient 
solution. Then 


(a) x* is optimal to BCMP if and only if X 2t = T 2 „ at x*. 

(b) At x if X 2 | < T 2 ) then f\ < /*; similarly, if X 2i > T 2I then /, > /*. 


PROOF: Part (a) Sufficiency: At x* let X 2 | = T 2 \. By definition of r 21 , at x* the DM is 


indifferent to unit gain in ,/j and A f 2 units of loss in f 2 where A f 2 = 
other words 


bU 

a./i 


dU 

df 2 


In 


(18) [/,(x*),/ 2 (x*)l - [/,(x») + 8|,/ 2 (x*) - A J 2 ■ 8,] 

where denotes the indifference of the DM and 8) is the differential change in /j. When 
A 2 i = T 2 1 , in the neighborhood of x*, a differential change 8, in ,/j is accompanied by a 
differential change of -(A f 2 ) • 8, in ,/j, as per the definition of A 2) . Hence, points in the 
neighborhood of x* are equally preferred to x*; hence, x* is the local optimum. By quasicon¬ 
cavity of U, x* must be the global optimum too. 


Necessity: Let x* be optimal to BCMP but X 2 ; 5* T 2] say X 2 , < T 2X . Around a neighbor¬ 
hood of x*, 3 a point l./j(x*) + 8 t , f 2 (x') — 8iX 2 )] where 8, is a differential change in ,/j as 
per the definition of X 21 . At X 2t < ^i, 

(19) (/ 2 (x*) - 81 X 2 «) > (/ 2 (x*) - 8,r 2l ) for 8) > 0. 

Since U is strictly increasing in both f\ and / 2 , inequalities (18) and (19) taken together imply 
that 

(20) l/,(x*> + 8|,/ 2 (x*) - 8|X 2 |] > |/,(x*) + 8|,/ 2 (x*) - 8| 7~ 2 ,] 

— [/|<X*)./ 2 (x‘)l. 

By transitivity of preferences 

(21) t/,(x*) + 8,. / 2 <x*> - 8|X 2 |1 > l/,(x*)./ 2 (x*)]. 

Obviously, (21) contradicts the optimality of x*. Hence, X 2! < T 2[ . Similarly we can show that 
x 2 | > r 2i . Hence the proof. 


Part (b) Let X 2 | < T 2{ . By (21) an increase in f x from its current value leads to a pre¬ 
ferred solution. Hence, an increase in ,/j is a locally increasing direction for U along the 
efficient frontier, Also x being an efficient solution is a solution to program (P v ) by Theorem 2 
and f\ (x) - v by Theorem 1. We noted earlier that x is a solution to program (P3) also, i.e.. 
Max U(\,f 2 ). Theorem 5 established the unimodality of g(v) which is the "value" of 
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optimal solution to program (P3); as /j(x) - v, the value of optimal solution to program (P3) 
is unimodal in f\(x). By the property of unimodality, U increases along a locally improving 
direction continuously up to the optimal solution and continuously decreases in the other direc¬ 
tion. Hence,/|(x*) > /i(x)or/f> /,. Similarly, we can prove the other case also. 

Note that Theorem 6 provides a convenient way of ruling out a portion of the efficient 
frontier. All that is needed is an interaction with the DM whether or not A 2 i - T 2 >. Depending 

> 

on his response we can narrow down the search of efficient points to only those with /, value 
greater than or less than the current /, value. We wish to point out here that the Surrogate 
Worth Tradeoff method due to Haimes, Hall and Freedman [9] is essentially based on a result 
similar to part (a) of Theorem 6. However, they do not have a result comparable to part (b) of 
the Theorem. Thus no "region elimination" is possible in their method. Also the DM has to 
specify how far T 2 \ is greater than or less than X 2) on a subjective scale of +10 to -10. 


We now formally state the procedure as follows: 

Step 1: Solve (PI) and (P2) and determine v* and w*. 

Step 2: Solve ((?*•) and determine /i(x*). Set v, /j(x 2 ). 

Step 3: Solve P y for v = v^ where v, < v^ < v u . Determine X 2 i, i.e., the Lagrange mul¬ 
tiplier corresponding to the constraint f\(x) > in the program P y . 

Step 4: Interact with the DM and present him with X 2 i and seek his comparative trade-off 
T-i i, i.e., whether T 2] S X 2! . Depending on the outcome of interaction eliminate a 

portion of the efficient frontier using the following rule: 

If X 2( < T 2 1 , update v, = v,,; go to step 5. 

If X 2 | > r 2 |, update v u = v^; go to step 5. 

If \ 2 | - T 2 \, stop; the current solution to P y is optimal to BCMP. 

Step 5: If |v, - vj < £ (chosen small number) stop; else go to step 3. 

Note that selection of v^ values can be done efficiently using the midpoint of the interval. 
For example, with 5 stages of interaction, we will be able to bracket the optimal solution to less 
than 4% of the original interval. In 7 stages, we will be within 1%, etc. Thus we have a finite, 
rapidly converging procedure using only "comparative tradeoffs." We shall illustrate the method 
using the sdme example given in Section 3. 

Numerical Example 

Solution Using Comparative Trade-Off Method: 

Step 1: As before v* - 60; w* - 3.2. Hence v„ - 60. 

Step 2: As before v, - 8. 
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Iteration 1 

Step 3: Using the 'Bisection' method (i.e., midpoint values) to choose v A values, v,, = 34. 
x A which solves P v for v - \ A has / values of (34,2.057); also the Lagrange multi¬ 
plier A 21 = 0.0571. 

Step 4. Interact with the DM and seek his comparative estimate of local trade-off r 2 |, i.e., 
whether T 2I > A 2) , T 2l < A 2 | or T 2l — A 2 |. Assuming that the DM is guided by 
the implicit preference function 


U(f,J 2 ) = fl n f 2 \ we get 

= 2/3 A," 1 ' 3 f 2 and ~r = f} n . Evaluating at x A , 
0 /| 0/2 

~ “ ~ (34) _ 1/3 (2.057) = 0.4233 and 
of [ 3 


~ = 34 2/3 - 10.4951. Hence, T 2l 
Of i 


6 U I d U 

9 /, / 9/ 2 


0.4233 

10.4951 


- = 0.0403. 


Thus, A 2 i “ 0.0571 > T 2 \ = 0.0403. Hence, update v„ to = 34 

Step 5: As the limits v, and v u are not close enough, the search has to be continued 
further on the reduced efficient frontier. The details of the calculation for further 
stages are summarized in Table 2. 


TABLE 2 — Solution of Bicriteria Problem Using 
Comparative Tradeoff Method 


Iteration # 

“i 


— 

V 4 1 

1 

(v A .g(v A )) 

*21 

r 2 , 

Remain Interval j 

(% of original interval) i 

1 0 

8 

60 


- 

- 


100 

1 1 

« 

34 

34 

(34,2057) 

0.0571 

00403 

50 

2 

2' 

34 

21 

(21.2.800) 

0.0571 

00889 

25 

1 3 

27.5 

34 

275 l 

(27 5.2.429) 

0.0571 

0.0589 

12.5 

1 4 

27.5 

30.75 

30.75 j 

(30 74.2.243) 

0.0571 

0.0486 

6 25 | 


5. BICRITERIA OPTIMIZATION APPLIED TO THE U.S. AIR FORCE 
CARDIOVASCULAR DISEASE CONTROL PROGRAM 

In this section we will demonstrate an application of the paired comparison method to a 
bicriteria optimization problem of designing a cardiovascular disease control program for the 
U.S. Air Force personnel, this part of a study currently undertaken by Purdue's School of 
Industrial Engineering, under contract with the USAF School of Aerospace Medicine. 

The United States Air Force (USAF) is planning a comprehensive program for reducing 
the incidence of cardiovascular disease (CVD) among its active duty personnel. This program, 
known as the Health Evaluation and Risk Tabulation (HEART) program, would consist of 
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(i) an extensive educational program about CVD, its causes and prevention, 

(ii) a risk screening program to identify those AF personnel who have a high susceptibil¬ 
ity for developing CVD, and 

(iii) a risk intervention program to treat such high risk personnel and reduce their CVD 
risk by medications and/or behavior modifications. 

In order for such a program to be viable, the expected benefits resulting from reduced 
incidence of heart disease must outweigh the cost of the HEART program. 

Risk Identification: A general cardiovascular risk profile, based on the Framingham Study 
[111, is presently being considered as the means of identifying high-risk personnel for treat¬ 
ment. With this approach a logistic regression model would be applied to each individual and 
would yield his CVD risk; namely, the probability of manifestation of cardiovascular disease in 
the individual in an eight-year time period. 


Risk Intervention: By rank ordering the CVD risk scores of USAF personnel, a certain 
percentage of the highest risk individuals are selected for specialized face to face therapy for 
modifying their elevated values of risk factors. At present, those high-risk individuals with 
elevated systolic blood pressures will be treated with medication, while those with elevated 
cholesterol levels will be treated with diet modification only. Both individual and group thera¬ 
pies will be used to persuade smokers to quit smoking. 


The recurring cost of the HEART program, excluding the initial start-up costs, is 
estimated between 7 and 9 million dollars a year depending on the total size of the therapy 
group. Hence, before implementing the HEART program on all airforce bases, the Air Force 
was interested in answers to the following question: 


Given the total size of the therapy group, determine the optimal 
number of Air Force personnel to be selected in each age group 
of flyers and nonflyers such that the total "cost-effectiveness" of 
the HEART program is maximized. 


The therapy selection problem has indeed two competing objectives. One objective is to 
minimize the USAF CVD risk and the other objective is to minimize the HEART program 
cost. By increasing the size of the therapy program, the CVD incidence in the Air Force will be 
progressively reduced but this can be achieved only at the cost of increased HEART program 
expenditure. If we follow a policy of selecting the high risk personnel from each age group for 
risk intervention, the graph of CVD incidence versus therapy size will be a monotonically 
decreasing function. Hence, there is an optimal expenditure level for the HEART program, an 
optimal level of CVD incidence, an optimal therapy size and an optimal selection strategy for 
that therapy size. Determination of all these quantities is a nontrivial problem for which we 
developed the following bicriteria mathematical programming model which was solved using the 
paired comparison method developed in Section 3. 
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BCMP Model 

Minimize 

Minimize 

Subject to 

Where: 


f\(x) 

fl(x) 

17 

I x, < S 

7-1 

a t < x, < b, 
j - I, .... 17. 


x = selection strategy (x|,x 2 .x 17 ) for each age group of flyers and nonflyers. 

/i(x) —USAF CVD risk for a given selection strategy x 
/ 2 (x) “HEART program cost for a given selection strategy x 
S — Maximum size of the therapy group (fixed) 

a f — Minimum number to be selected from age group j (flyers and nonflyers); 

j = 1.8 correspond to flyers in age groups 20-24.55-59, respecitvely; 

7 = 11.17 correspond to nonflyers in age groups 15-19.55-59, respectively. 

b, - Maximum number that can be selected from age group j. 

Even though / 2 (x) could be determined analytically, f t (x) is indeed a random variable 
because of the unpredictability of the effectiveness of the risk intervention programs on modi¬ 
fying the risk factor levels. Hence, a comprehensive simulation model known as P-HEART 
(Purdue Health Evaluation and Risk Tabulation) was developed by the authors [14]. P- 
HEART simulates the USAF population and performs risk identification, therapy selection and 
risk intervention. The description of the simulation model, its validation and findings are given 
in Reference [14]. Through extensive simulation runs, the expected value of ,/j(x) for different 
selection strategies were estimated and a polynomial equation was fitted to get analytical forms 
for f\(x) by age group and flying status 

The paired comparison method was then programmed within the frame work of interac¬ 
tive solution of BCMP for minimizing CVD risk and HEART program cost. The Generalized 
Reduced Gradient Algorithm was used to solve the nonlinear programs in Steps 1, 2, and 3 of 
the interactive algorithm. The interactive program was extensively tested at Purdue using Air 
Force Officers as decision makers. The program has been well received by the policy makers in 
the USAF HEART program office and is being used as part of the long term project planning 
currently underway in the Air Force. For more details on the interactive computer program 
and its output, the reader is referred to Sadagopan [16]. 

6. CONCLUSIONS 


In this paper we have developed interactive procedures for the solution of bicriteria pro¬ 
grams. Even though the entire set of efficient solutions can be generated and pictured, such 
iteractive procedures do have merit due to the following reasons: in many real life problems 
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complete generation of all efficient solutions may be computationally expensive; even granting 
such generation, the DM is unaided in the selection process. Hence, the procedures developed 
in this paper may be thought of as "decision aids" [15], 

The primary power behind our procedures is the unimodality property (Theorem 3). A 
similar result is available in Geoffrion [7], However, our proof of unimodality is entirely new 
and allows extension to more general cases unlike Geoffrion's result. Also our relationship 
between Lagrange Multipliers and local tradeoffs is new. For more details on the bicriteria 
results, U.S. Air Force application and extensions to the general multi-criteria problem, the 
reader is referred to Sadagopan [16]. 
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ABSTRACT 

There are a great number of queueing xyxtemv including the A/' A/’/c. the 
G/'/A//c ami the discrete Ch in I queue in which th - te probabilities are 
determined h\ repeated queue equations This rape- . a simple, efficient 
and numerical!! stable algorithm to calculate the sta i >b. .ties and measure 
of performance for such systems The method avoi th complex artth- 
metric and mainx manipulations 


1. INTRODUCTION 

This paper gives a new method to find the equilibrium probabilities of continuous-time 
Markov processes in which the rale of going from state in to n is equal to J„-„, for any n 
exceeding a certain limit, say r. The largest and smallest possible jumps are It and —g respec¬ 
tively. that is, the subscript of d k runs from -g to h. These concepts will be further expanded 
in the next section. Here, we demonstrate it shortly, using the \i x lKVi 1 queue as an example. 
In this queue, arrivals of size /, 1 ^ ./ ^ li occur at rate K,. and departures of groups of size i 
occur at a rate of /x,. 1 ^ < g. Hence. 

d, = A,. I <./</•; 

z/_, = fi,. i ^ n, 1 ^ ^ g. 

The paper also considers discrete Markov processes with a similar structure. In this case, 
there is a probability of p k to go from state n to n + k. Queues with such structure include the 
system size of the GI*/M/c queue (41 and the waiting lime of the discrete G7/C/1 queue (5j. 

At present, there are several approaches available to solve such queueing problems. In 
the first approach, one sets up generating functions, and uses partial fraction expansions to 
invert these generating functions. These expansions require one to find the roots of an equa¬ 
tion of degree g + h. The same is true if one uses operators on the queueing equations and 
determines the characteristic roots. I or practical applications, these approaches have some 
disadvantages. The equations are often difficult to solve, especially when they are of a high 
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degree. Their roots are usually complex and difficult to work with. Moreover, there may be 
double roots which lead to further complications. Neuis (3! therefore suggested solving these 
problems using matrix equations. 

The method in this paper combines the approach of N'euts wiih the classical method. In a 
sense, it can be considered as a modification of the method of Neuts. The main difference 
between our method and the Method of Neuts is the following. In Neuts - approach, a matrix R 
of dimension h x h has to be calculated. In our approach, only the first column of R is calcu¬ 
lated. which reduces the number of operations approximately by a factor of h. Even though the 
method is similar to the one of Neuts, its derivation is done, using the classical approach. In 
this way, our method bridges the gap between the classical approach and the approach of Neuts. 


2. PROBLEM DEFINITION 


In this section, a precise problem definition is given and demonstrated, using a number of 
examples, first, we discuss the discrete case. Let there be a Markov chain with the transition 
probabilities. tun ^ 0 The equilibrium probabilities of this chain are denoted by tt„. 
that is, n„ is probability to be in state n As is well known, the tt„ are determined by the fol¬ 
lowing equilibrium equations, provided the system has an equilibrium 

(I) IT,, = £ it m n ^ 0. 

We assume that there is a certain limit /-such that for n > r, p m „ = only depends on the 
difference between m and n In other words, if Q represents the state. d k is the probability to 
increase Q by k from n to n + k. If tr„, is defined to be zero for m ^ 0. equation 11 ) becomes 
therefore for n ^ r 


tt„ =* £ TT m d n _ m = £ 7r„_.*<4, n 2* r. 

m —oo A—°o 

We furthermore assume that d k = 0 for k < -gor k > h. Thus, one has 

h 

n " ~ X W 'i-* d k‘ " > r. 

*—x 


If one solves this equation for n one gets 

h 

(2) n„ = £ n " a e k- " ^ r - 


The e k in this expression are defined as 


(3) 


e k = d k/l I ~ d 0 ), k t* 0| 
e„ = 0. J 


for n < r. one uses equation (1), that is 


14) 7r„ — £ n m p m „. n < r. 

m-n 

Equation (2) will be referred to as the queueing equation , whereas equation (4) gives the initial 
conditions. There are r initial conditions. However, it can be shown that one of them is redun¬ 
dant. In its place, one has the normalizing condition 
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<s) 

n-n 

If there exists a unique equilibrium vector ( 77 ,,!. this vector is fully determined by the queueing 
equaiion. the r — 1 initial condition and the normalizing condition. 01 course, if r — 1. there 
are no initial conditions. As is well known, an equilibrium solution will exist if 

p = X / <tj x, i ii , < 1 

This paper will show how to find this equilibrium solution in an efficient way. 

In continuous-time Markov chains, the equilibrium probabilities are uciermined by 
0 = X v,„„. n ^ 0. 

Here. q m „. m ^ n. are the transition rates, and 

dm ,,, = ~ X ‘i’" "' 

We now consider problems where 

dm,, = n ^ r. 

Here, </„ is defined as 

= - X 4 • 

As in the discrete case, it is assumed that <l k = 0 outside the range —g ^ k ^ h. If 77 ,, = 0 for 
n < 0, one finds thus from equation (6) 

h 

0= X Ttn-kdk- n > r. 

k — K 

This equation can again be solved for 77,,, which gives 
/; 

(7) 77 ,, = X n„^ k e k . n ^ r. 

k — K 

Here, 

e„ = 0 

<8) t\ = d k l (—(/,,), k * 0 . 

Note that equation (2) and equation (7) have exactly the same structure. In particular, > 0 
and the sum of the e k is one in either case. The normalizing condition is given b\ (5) both in 
the discrete and the continuous case. The initial conditions of the continuous case are easily 
found from (6) by letting 0 < n < r. 

We now give 3 examples which will be used later for the purpose of demonstration. 

EXAMPLE 1: First, consider the waiting time in queue of a discrete G//G/1 queue. Fol¬ 
lowing Ponstein (51, a Markov chain is set up, in which the state space are the waiting times 
rather than the number of elements in the system. To set up this chain, it is assumed that the 
difference between the service time (S) and the interarrival time (A) is (a) always integer and 
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(b> always between —.y and h. It is well known (see e.g. Ponstein (51). that the waiting time of 
customer / exceed the one of customer / — I by an amount S — .4, provided, of course, that 
customer / - I had a waiting time above — (.S' — 4 I. Consequently, if H represents the wait¬ 
ing time of a customer in equilibrium. W increases by A with probability P(S — .4 = A). This 
means 

</* = P<S - A = A), < A ^ /; 

and 

e k = P(S - A = A)/[I - P(S = .4 )|. A * ()) 


In the problem considered, there is only one initial condition, and this condition can be 
dropped as redundant. The queueing equation, together with normalizing condition is thus 
sufficient to determine all tt Here, tt„ represents the probability of having a waiting time of 
n. 


EXAMPLE 2: In the Gl h /M/\ queue, arrivals occur in groups of size h. Let a, be the 
probability of servicing / elements between two consecutive arrivals and let n„ the probability 
of having n elements immediately before an arrival. It is easily verified that the queueing equa¬ 
tion is given by (2) with r = I and 


( 10 ) 


(\ = «/,—*/< 1 - «„>. -oo < A < /;. A * 0] 

<■'<) = j 


Here. is the probability of serving A customers between two successive arrivals. The initial 
conditions consist only of one equation which can be omitted. Also, in the case considered, g 
is infinite. However, for practical calculations, it is sufficient to carry only a finite number of 
<V 


EXAMPLE 3: Consider the queue. A ( , I < / < h. is the rate at which groups 

of size ./arrive. No service is done unless at least g elements are in the system. If the number 
in the system is g or more, there is a rate of /*,- 1 ^ < g, that / elements leave the system 

after having received service together. For this problem, one has the following transition rates 

q„, „ = A„_ m , 0 < n — m ^ h 

q m „ = n m -„. 0 < m — n ^ g, m > g. 


All other q m n . m ^ n are zero, in particular, all m < g , n < m are zero because nobody 
leaves if there are less than /.'elements in the system. The q mm become 


q mm = — £ A, m < q 

/-i 
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A little reflection shows that r = and the queueing equation is given by (5) with 



f " * ; 


\ - A A. / 

Z x ,+ 2 > 

h K 

, k > 0 

’k = Ma/ 

Z + I>> 

/-I /-I 

, A < 0 


The initial conditions are given by (7), that is 

oo h 

0 = £ = -Vn £ A, + 

it-0 /-l 

0 = £ = Z ”V-/ X / ~ Z X ' + Z « < £• 

a- ft i-1 /-i (-i--/) 


Again, the convention was used that tt m = 0 if w < 0. From the initial conditions, one can be 
omitted, and it turns out to be convenient to omit the one corresponding to n = g — 1. 


Above, we discussed a few examples which are amenable to our method. Additional 
examples can easily be generated. We now show how to solve these examples. 


3. THE MAIN RESULT 

The fundamental result of this paper is the following. If p < I. and if the Markov chain 
is evgodic. there are h values o, > 0, 1 < ./ ^ h such that 

h 

( 12 ) tt„ = £ a,TT n . r n > r. 


The a r 1 < j ^ h are determined in a unique way by the e ks ^ A < h. Further¬ 
more, 

h 

!«,<'■ 

./-i 

This paper will present several methods to find the a r Once the o, are found, the problem is 
actually solved. Together with the r — 1 initial conditions and the normalizing equation (12) 

gives exactly r + l equation 10 determine n n , . n r . Moreover, the generating function 

of the 7 t„ can be expressed in terms of o, and 7r () . 7T|.7r,_,. To see this it is convenient 

to define a n as -1, and set 

h 

A (z) = £ a z ' 

,-n 

P(:)= I 7r„r". 

In this case, one finds, as will be shown in Section 4: 

(13) A<z)P<z)- £*,-, 2 *. 

,-n n-, 
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When setting r = 1. Pi:) becomes I, and 

h r -1 f— I <■ 

.4(1)= £ a, £ n-„_ / = £ 7r„ £ a,, <' = mini/- — n — 1. //). 

/-ft n-i »-0 /—0 

This means 

/> r- I 

(14) £ a, = — 1 + 0 | + ai + ... + a,, = £ 7r„ £ a,, e = min(r — a — 1, /;). 

/-a n-n ,.n 

This equation provides a convenient normalizing condition and can be used in place of (5). To 
find L. the expected number in the system, one can differentiate (13) and find 

h r— I 

A'iz) Pi:) + .4 (r) P‘i:) = £ a, £ 

,-n „~i 

If r = 1, Pi:) = 1 and P'i:) = L Using these values, the above equation gives after some 
calculation 

(15) L = £ a, / £ a,. 

/-0 ,-o 

Equations (14) and (15) simplify if r = 1. Then 
/I (e) P(e) = ao 7r « = “ 7 r o- 
If one sets r = 1, this gives, 

(16) ir 0 = -Ai\) = 1 - £ a,. 

/-i 

This means that 

I a, < 1. 

If Q is a random variable representing the state of the system. EiQ) can be found as: 

(17) EiQ) = £ nir„ = -A'i\)/Ai\) = £ ./a,/ tt,,. 

«-0 /-I 

The ideas just presented shall now be demonstrated, using the examples discussed above. 
First, consider the discrete GI/G /1 queue. Specifically, suppose that the distribution of interar¬ 
rival time and the distribution of the service time are as given in Table 1. 

TABLE 1 — Distribution of the Interarrival Time 
and the Service Time. 


k 

1 

2 

3 

PIS = k) 

0.4 

0.3 

0.3 

PiA = k) 

0.1 

0.3 

0.6 


Using these values, one can find the e k , -g ^ k < /t, from equation (9), and the e h in turn, 
determine the a,, 1 < j ^ h (see Sections 4 and 5). For the problem described by Table 1. 
one finds, 

a, = 0.23156 
a 2 = 0.05039. 
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Equations (16) and (17) give now 

Ti-n = P( W = 0) = 1 — a, — a 2 = 0.71804 

E(W) = ( fll + 2 a 2 )/P( W = 0) = 0.46287. 

Since 77_, is assumed to be zero, one can use (12) to calculate tti , 7r 2 , w-, and so on. 
77 | = 7T()<? [ + 7T_|fli = 0.1663 
77 i = TT i ZZ i + 77|,a 2 = 0.0747 


Since all a, > 0, this recursion is numerically very stable. 

Next, consider the Gl 1 '/M /1 queue. Specifically, let the arrival time distribution be deter¬ 
ministic, let h = 2, A 2 = 0.4 and n = 1. The a k of equation (10) are in this case Poisson- 
distributed with parameter 2.5. One finds in this case 

a, = 0.41045 
a 2 = 0.13711. 

Again, all entities of interest can be calculated. L. the expected number in the system preceed- 
ing an arrival becomes for instance 

L - (a, + 2a 2 )/(l - a, - a,) = 1.5132. 


Finally, consider the M x /M y /\ queue of Example 3. Specifically let 
= \ 2 = Mi = M: = M.t = *• 

In this case, the a, turn out to be 
a, = 0.34960 
a 2 = 0.24582. 

The initial conditions are now 
0 = —2no + 77 1 
0 = 77 0 — 277 1 + 77 3 + 77 4 , 

773 and 77 4 are given by (12) as 

773 = a]77 2 + a 2 77 | = 0.34960t7 2 + 0.24582t 7 t 

774 = a|77 3 + a 2 77 2 = (a | 2 + a 2 )t7 2 + aia 2 77| — 0.3680477 2 + 0.0859477,. 
Using these values, the initial conditions become 

0 - 2t7 0 + 0.2458277, + 0.34960t7 2 

0 - 77 0 - 1.6682477, 4- 0.71764t7 2 . 

As a third equation, one uses the normalizing condition as given by (14). 

— 1 + a, + a 2 — i7n(~1 + a, + a 2 ) + tt, (—1 + a,) — 77 2 
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or 

0.40458 = 0.40458tto + 0.65040^, + ir 2 . 

These 3 equations ean be solved for ttd. it i and wd. 

The result is 

7r„= 0.06741. a-, = 0.15840. n 2 = 0.27428. 

It is now simple to find L from Equation (15). One has 

L = [— <i + 2 ttO + a | (ff,, 4- 2 tT| — 1) + a 2 <2-jt 0 — 2)]/(— 1 + a | + a 2 ) 
= 3.4128. 


Thus, once the a, are found, the problem is solved easily. The next two sections show 
several methods to find the a,. 


4. THE SIGNIFICANCE OF a, AND THEIR CALCULATION 


In this section, it is shown that equation (12) holds, and how one can find the a, in the 
general case. From the initial conditions and from the queueing equation, one finds the follow¬ 
ing expression: 

P(:) = U(:)/V(:). 

Here U(:) is a polynomial of degree r + r — 1 which is of no further interest, and l'(r) is 
equal to 

it+h 

VI:)- I fn:*-:' 

*-n 

It is well known (see e.g. 15]) that U(r) has a zero at : — 1. r — 1 zeros inside the unit circle 
and h zeros outside. It follows that K(e) can be written as 

Viz) = A (z) Biz). 


Here 


Aiz) ^ — 1+0|7 + a 2 z 2 + ... + a h z h 


contains all the h zeros outside the unit circle. 


Because of the convergence of Piz). Uiz) must be divisible by Biz). Since Uiz ) has a 
degree of r + r - I, and Biz) has a degree of r. Ciz > has a degree of r - 1. Consequently. 

r-\ 

Uiz)/Biz) = Ciz) = £ c,r'. 


or 

Piz) = Ciz)/A iz). 
This equation can be written as 
(18) Piz)A(z) - Ciz). 
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The coefficient of of the product fMrM (r) equals, as is well known 
£ °. 7I Vr 
for n ^ r - 1. this means 

(19) £ <i,7T„_, = C„. 

/-II 

For n > one finds 

(20) £a,ir„_, = 0. 

,-o 

Equation (20) is identical with Equation (12), and Equation (18) and (19) can be combined to 
give Equation (13). This proves Equations (12) and (13). 

From the above discussion, the significance of the a , becomes clear, TO) can be factor¬ 
ized into two factors ••)(;) and Biz), and the a, are merely the coefficients of A <r>. provided 
£j () is defined as -I. The problem is thus reduced to finding A(z). The most convenient 
method to do this seems to be the following one. One starts with the identity, which will be 
proven in Section 6. 

a, = e, + £ e_, «, ,, 1 < j < h. 

t-1 

at),i = a, , I <./'</» 

a,+i., = a, i a, + a, l+[ . I < ; < jtf - 1. I < ./ < /; 

Furthermore. a u is to be taken as zero fo. j > h. Equations (21) and (22) can be used to find 

a, by successive approximation. As starting values, one can use a," - e r j = 1.2. h. 

These values can be used to replace the a , on the right of (22). and a new approximation for a, 
can be found from Equation (21). In this way, one continues uniil a suitable Slopping criterion 
is satisfied. 

We tried two different stopping criteria. First, we required that the n, generated by Equa¬ 
tion 02) satisfy the queueing equation with a specified precision e. It was found that for high 
values of h and/or #, this criterion performed poorly. The reason for this is simple. If h is say 
100. a small change in t/| IK > will have a dramatic impact on p and, consequently, on the perfor¬ 
mance measures of the system. As an alternate stopping criterion, we used the change of the 
calculated /f'(l) = I ja / to decide when the precision of the a, is adequate. Since AA\) is 
closely related to the expected number in the system (see Equation (17)) this seems to be a 
good stopping criterion. 

If there are q iterations, one needs 4qgh operations to find all a r The reader may want to 
verify this. If one adds all e_,a,, to the sum at the right of Equation (21) as soon as the a,, 
becomes available one only needs to store the a,, for the current value of i. If this is done, the 
algorithm requires only an array area of 4 h + jzas the reader may verify. The algorithm is thus 
efficient and does not require huge arrays. 

The algorithm described above was programmed in FORTRAN and run on the DEC 
2060. The execution times were negligible. A problem with h = 100 and g - 100 and p = 0.9 


( 21 ) 

with 

( 22 ) 
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<>nlv required 8 4 seconds and 53 iterations until all a, were obtained. At this point there was a 
relative change of 4'< I) of less than 0.0001 between two iterations. The computer time to find 
i be u lor e = 100 and h = 5 is 0.2307 seconds (19 iterations). The method performs thus 
vctv well 

V Al.TKRNATF. METHODS 

In gain a check for our results, and also to get further insight into the problem, we 
decided lo trv other methods as well. These methods are based on l'(z) = V(z)/(z — 1), that 
,‘is I i.-i was deflated by z — I l'(c) can be written as 

1 <;> = 4 (z)fl(z). 

where /f<:) is equal to fi<z)/(z - 1) This deflation decreases the degree of the polynomials 
one works with by one. and it also increases the difference between the zeros of A (z) and 
H i.-). improving thus the efficiency of the algorithms to be discussed. A further deflation by 
the onlv positive /ero of .4 (z) is possible. 

To factorize 1(c). two methods were used. The first one is given in [ 1, page 1581. The 
second one used the fact that f'(z) can be interpreted as the characteristic polynomial of a 
difference equation If these difference equations are used to calculate values x„, x.+\, x „+2 
recursively, ihe zeros greater than one will dominate, and eventually, the effect of B(z) will 
become negligible The ,v„ are thus almost identical to a series generated by a difference equa¬ 
tion that has 4(c) as characteristic polynomial. In other words, the ,v„ will almost satisfy the 
following relationship 

h 

(23) ,v„ = £ a, x„+ r 

If x„ is known for n = A, k + 1 . k + 2h - 1, this gives h equations for the h unknown 

a r 


The ,v„ will satisfy Equation (23) precisely if the initial conditions x n , X\ .satisfy 

(23). Consequently, one can repeal the above algorithm several times lo gain a higher preci¬ 
sion for the a r In each repetition, one uses the a, obtained in the previous iteration in order 
to calculate the starting values needed. In the first iteration, one can use a , n — d r This algo¬ 
rithm gave good results for most cases we tried 

The disadvantage of the two algorithms just mentioned is that they require the solution of 
h equations in h unknowns during each iteration, and this requires h 3 operations. This means 
that they are inefficient as compared to the algorithm suggested in the previous section, at least 
if h is high. 

Of course, A (r) can also be found, provided one knows all the zeros of V(z) outside the 

unit circle. Let Z].Z 2 .z 3 . z h be these zeros- A (z) becomes, given one uses the fact that 

a o - -1 

A (z) = — (Z — Z\) (Z — Z 2 ) . . . (Z — Z h )/U— Z|) (— z 2 ) ... — ( — z h )\ 

= — (1 — z/Z|) (1 - z/z 2 ) ... (1 - z/z h ). 

The factors can now be multiplied in the usual way, (see e.g. [4]) giving 
A (z) — — 1 + a\ z + a-iz 1 + ... + a h z h . 
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Note that this procedure even works for multiple roots, provided, of course, all are known. 

A final comment shall now be made. Suppose one has the h + k probabilities ir A + |. 

n k+: . n k Then, one can theoretically calculate nm > k + h + it', recursively by 

using Equation (2), that is 



This recursion is numerically unstable, that is, the round-off errors will increase with each itera¬ 
tion. The characteristic equation of this difference equation is Ml/.v) = ,4(l/.v) l/.v). and 

any recursion based on it will eventually be dominated by Z?(l/.v). (Indeed, this very effect 
was used earlier to find A(:) s Instead of doing such a recursion, it is better to use 2h subse¬ 
quent it,,, and use (12) to obtain It equations for the a r j = 1.2. h. 

6. THE SIGNIFICANCE OF THE a,, 

This section proves Equations (21) and (22) and establishes the relationship of our 
method with the method of Neuts. We start to prove the following equation, in which the a ,, 
are calculated as given by (22) 

h 

(24) n„ +l = £ a, , jt„_, • / ^ 0. n > r 

./-i 

Since a,), = a,, this equation is certainly correct for / = 0. Moreover, if it is true for /, it is 
also true for / + 1. One has, if a hh+l = 0 

/,+1 h* I 

”n+,±i = J) a u n » + 1-, = «,.t + £ a ‘-' 7r » + 1-/ 

1 i-2 

It It 

= a, i £ a,na, . ,+| tt„_ / 

/- T i-i 

it 

= £ [a, | a, + a,n. l ]n„- r 

i- 1 

Because of (22), this gives 

h 

7T n*t*\~ £ a i+\.i * n-r 
I- I 

This proves that Equation (24) is correct for / + I, given it holds for i. To prove Equation 
(21). one rewrites the queueing Equation (2) as 
i i ' 

tt* - £ e.,n„ +l + £ e,TT 
(-1 /-I 

Using Equations (12) and (24) this gives 

h h x h 

L - £ <’,"•«-/ + Z t ’-' Z 

,-t ,-t ,-t i-\ 

= L ^ + L c ’-< a J w «-r 

,-l ,-t 
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Consequently, 

a, = e, + £ e_,a,,. 

/-i 

and Equation (21) is proven. 

Equation (24) can also be used to simplify the initial conditions. Consider for instance 
Equation (4) which can be written as 

°o r— I 

77 „ = £ TtmPm.n = ^ 77+ £ Pr+m.n™ r+m- " < Z. 
m-() m-0 m-l> 

If 77, +m is replaced by (24), one obtains after some calculation 

r-/,-l r— 1 oo 

77„ = 77 + £ 77<7* ,_ m l, 0 Sj H Sj /'. 

m-n m-r-l, A-ll 

Together with Equation (14), this provides a set of a independent equations in the r unknowns 

TTo. 77 i.77,— |. 

We now compare our results with the ones obtained by Neuts. [3,41 Neuts sets up the 
following matrix equation for the unknown matrix R. 

(25) R = £ R"A„. 

«-o 

The ,4„ are problem-dependent matrices. In our case, they become 

A„ = [a;',\. « - 0.1.2. Ig/h 1 + i 

Here, [#/);l is the lowest integer above n/h, and 

•4,”, = G, 

All other /!„ are zero. As before, all e k outside the range -g < A- < h are assumed to be zero. 
If one defines 

~ [77|. 77*/, + 2 . 77*/, + /,I, 

one has according to Neuts 
77* + l = TT k R. 

If R = [/•,,), this gives 

h 

ff (*+t>i)+/ *■ Y* n *fr+/ r /.< 

/-i 


"^n + i = I- i r li+\-i.i • 

/-I 

When this equation is compared with (24), one finds 

a u = 0 / + I-/./ + I • 

In particular, 

a o.i ” Oj ” Ot+1—/. i * 
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The a t are thus eoniained in the first column of R We also note that liquation (21) follows 
from Equation (25). as the reader may verify. 

To calculate R according to (25) one has to find the terms R"A„. n = 2.3 . ln/h\ + I 

and this can be done in 2 (#//>] matrix multiplications. However, each matrix multiplication 
requires 2 h‘ operations, which means that there are 4 h*lg/hi. or approximately 4 /re opera¬ 
tions per iteration to find R. not counting the matrix addition Thus, by taking advantage of 
the special structure of (he problem one obtains a more efficient algorithm, namely the one 
described in this paper. 

The equivalence of the a, with r h + 1 _,| is of great importance. In particular. Neuts (21 
proved that all r,, are nonnegative. Consequently, all a, are nonnegative Neuts also finds that 
the r, ! have a probabilistic interpretation. Indeed, the r, , are the expected number of visits to 

state n + H + i starting at n + I without hitting In + I. n + i . n + /;). This means 

that a, gives the expected number of visits to slate n + / starting at n without hitting 

I«+./—/».«+./'- 11. 

7. CONCLUSIONS 

This paper derives an extremely efficient algorithm which can be used to numerically 
determine the slate probabilities of many queueing problems efficiently and precisely This 
algorithm was obtained a by combining the classical approach with the approach of Neuts. 
Indeed it bridges the gap between these two approaches and opens new horizons for further 
research. 
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INDIVIDUAL HEADSTART STRATEGIES 
FOR COMBATING CONGESTION 
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ABSTRACT 

This paper examines the process by which a user of a queueing system 
selects his arrival time to the system to compensate for unpredictable delays in 
the system if he wishes to complete service at a particular time Considering 
the case in which all the system users have already decided on their arrival 
limes to the system and will not change these times, this paper investigates how 
a new user of this system develops his strategy for selecting his arrival time 
The distribution of this customer's arrival time is then obtained for a special 
case. 


1. INTRODUCTION 

If a user of a queueing system can choose when to enter the system, then he would prob¬ 
ably choose to enter when the delays he experiences will be minimum. However, if he wishes 
to complete receiving service at a particular time, called the "target time," then he would have 
to arrive at the system early enough to allow for his service time and unpredictable delays in 
the system. This allowance was termed "headstart" by Gaver [31 and the same terminology will 
be adopted here. The strategy by which a user selects his headstart is studied in this paper. 

Consider the case in which all the system users already have fixed headstart strategies and 
will not change these strategies under any condition. This paper investigates the case in which 
a new user who wishes to use the system develops a headstart strategy. This will be done by 
assuming that this user, whom we name U (cf. Alfa and Minh [2]), attaches some perceived 
costs to delays and to early and late completion of service, and that he selects his headstart in 
order to minimize his total cost. Gaver [31 obtained the "best" headstart strategy for this user, 
in this type of situation, on the assumption that the user wishes to minimize his expected total 
cost. This paper proposes the headstart strategy usually adopted by such a user, on the assump¬ 
tion that he selects his headstart each day in order to minimize the total cost he incurs. It is 
assumed that on the first day he uses the system he does not know what the delay distribution 
is and therefore selects his headstart arbitrarily, and hence probably incurs high total cost. 
However, for the following day, he uses his knowledge of the previous day's outcome to try 
and reduce his total cost; and this continues until a stage is reached in the long-run when he 
settles for a particular headstart and does not change any further. Let us call this stage 
"steady-state." This headstart reflects his arrival time at the system and our interest is in the 
distribution of his arrival time at the system at steady-state. 
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This paper does noi restrict its analysis to nonrush hour traffic or linear cost parameters as 
done by Gaver (3], 

2. PROBLEM FORMULATION 

2.1 Assumptions and Definitions 

Consider a queueing system used by a finite number 1,(1 > 0), of customers who behave 
identically and independently. For each day. let us observe the system at equally spaced epochs 
sequentially numbered 0,1,2,3.... and assume that all the customers arrive at and depart from 
the system only at instants immediately prior to these epochs. 

It will be assumed that once a customer is served he does not return to the pool of poten¬ 
tial customers. This assumption and that about the number of customers being finite are not 
necessary for the analysis but they allow us to use the results from Alfa [1] and Minh [4). 

Assume that each customer, considered separately, has the time-dependent probability 
of arriving at the system prior to epoch n + I, given that he has not arrived at the system by 

epoch n 0/ = 0,1,2.3_). Assume that the service times of the successive customers are 

independent, identically distributed random variable S in the set {1,2 . A/): M < °o. 

Suppose the new customer. U who also wishes to use this system, chooses to arrive at the 
system at epoch a = n. Let the delay he experiences by so doing be W". The distribution of 
W" can be obtained as shown in Alfa (!]. Note that the addition of this customer U to the sys¬ 
tem has raised the total number of customers to / + I. Although the distribution of i/'s arrival 
time is not necessarily the same as that of the other / customers, the results in Alfa [1] can still 
be applied keeping in mind that W" is a conditional random variable. The distribution of If s 
arrival time shall be developed in this paper. 

Let V" be the time spent in the system bv U. given that he entered the system at epoch n, 
then V" = W" + S. Let F," A Pr{V" = /). 

For practical purposes it will be assumed that U shall arrive at the system only at an epoch 
between 1 and ;V; where 1 < \ < 

Suppose U wishes to complete service at epoch r, where for convenience 1 < r $ ,V. If 
he arrived at the system at epoch a = //, got delayed W" = i units of time and took S = / units 
of time for service then he would complete service at epoch n + / + j. One of his objectives is 
for n 4- / + ./ to equal t. However, if /? + / + /< r or » + /+./> r then he is early by an 
amount r - n — i - / or he is late by an amount n + i + j - t. respectively. He attaches a 
cost to either of the outcomes and these costs are not necessarily the same. In addition, he 
attaches a cost to the time spent in the system i + ./. He would, therefore, incur a total cost 
C(n, i + j), where CO/, / 4-./) is the sum of the cost of time spent in the system and the cost 
for either earliness or lateness, depending on which is the outcome, keeping in mind that if 
« + /+./ = r then he only incurs the cost for his time in the system. 

Let C„ be the total cost incurred by U given that he arrived at the system at epoch n 

where C„ can only assume the values in the set {C(n.il|» — 1.2.(/+ II x M ). For 

brevity, C„ will be termed the cost associated with epoch n. 
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Let us define the following sets that will he used in the next subsection 


(it 

G A (i 1 i = 

1.2. ... 

. •/ + 1 t x A/) 

III) 

L A (/I/ = 

1.2.3. . 

.. . ,vl 

liiil 

L„ A (ili = 

1.2.3. . 

... n — 1. n + 

<iv ) 

L", A /.„ n 

/.„. 



If there are v epochs. (I < v ^ ,V - 2), LV > 2). labelled m t . nu .m, for which 

C m = C m| = C,„. = . . . = C,„ where Vm, € Z.",. r = 1,2.v. then let 

(v) /." , A (tit = ttt,. ttt>.in,, Vm, € L',l, . r = 1.2.v|. 

Let 

(vi) L n m , A — Z.," i.e. , is the set difference of L m and L ",,, 

2.2 The Model 

2.2 / The General Model 

Suppose on the i/” 1 day. (d = 1,2, ...I. U arrived at the system at epoch a‘'= n. Let 

a'J A Pr(a‘ y = t;(, tt € L. U would incur a total cost C„. and by definition 

(1) Pr(C„ = CU /» = V”. 

We assume that IT s choice of arrival time epoch for one day depends entirely on the out¬ 
come of the previous day's choice of arrival lime i.e., on the total cost incurred. V wishes to 

minimize his total cost. Hence, a 1 ** 1 is a decision variable with an action space L, and a J 
depends on the cost incurred on day d. Given the action on the <f h day the action for the 
(d + 1 ) ,h day, a‘ /+1 , is thus a Markov Chain. 

Let r„ m A Pr(a‘ y+I = m la 1 '= n) be the transition probabilities of the chain These transi¬ 
tion probabilities define IT s strategy. 

If s intention is to choose an arrival time epoch that minimizes his total cost each day he 
uses the system, hence for the (d + l) lh day, his strategy for choosing his arrival time epoch 
will be 

(2) t nn - Pr{C„ < C u ; V u € Zj; 

(3) i nm = Pr(C m < C B . Vm € L m \ 

+ V —~ Pr|C_, = C m < C„, Vm, € L" v ; Vm € L n m , h 
“T v + I 

Vm ^ n. 
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The arguments leading to (2> and (3) above are as follows: 

la) If L arrived at the system at epoch n on the if h day he will not change his arrival 
time to any other epoch on the id + !) lh day if he can not reduce his total cost by so 
doing. This leads to (2) 

(b) If i arrived at the system at epoch /ton the </"’ day he will change his arrival time to 
epoch w on the (d + 1 )' h day if this could reduce his total cost, prov ided there is no 
other epoch which has equal potential reduction in total cost as epoch w— this leads 
to the first term on the R.H.S of (3). In addition, however, if there are v other 
epochs, other than n and w. at which U could choose to arrive at the system and 
achieve the same reduction in to’.:.! cost as arriving at epoch m. then L would arrive 
at any of these v epochs or epoch w with the same probability, provided the total 
costs to be incurred at each of these epochs, including epoch w. is less than at any 
other epoch. This leads to the second term on the R.H.S. of (31 (2) and (3) above 

constitute the transition probabilities that define L~ s strategy 

Define an ,V x ,V matrix T = U„„l and an .V vector A ,/ = |of. a* . u‘vl- Then L 

selects his (</ + I )' h day's headstart according to 

(4) = A 1 ' x T. 

THEOREM 1: £/„,„=l. 

PROOE: Consider one epoch n. (n € L). and the associated cost C„. If we compare this 
cost C„ to the costs associated with other epochs in the set L„ then C„ is either the minimum 
cost, one of the equally minimum costs or neither of the two. This can be stated as 

(5) Pr{C„ < C u ; V u € L„\ + Pr{C„ <C„: V M €/.„}-!. 

Substituting i„„ for the first term on the L.H.S. of (5) gives 

(6) !„„ + Pr{C„ <C„; Vw €/.„)=!. 

The second term on the L.H.S. of (6) is the probability of having at least one epoch, in 
the set L„, whose associated cost is less than C„. This implies that there is at least one epoch, 
other than n, whose associated cost is minimum. If there is only one such epoch then the pro¬ 
bability of such event is given by 

£ Pr(C m < C„ :Vm- € L m \. 

y mfZ, 

If. however, there are v other epochs (I ^ v ^ ,V — 2). labelled m t . m : . m,. Vw, 6 L*. 

I < r < v, such that C,„ = C m = C,„, = ... = C w . then the probability of such event is given 
by 

v-z I _ 

£ £ —TT Pr < c m = C m < C,; Vw, € L", s . Vw € L'l v ). 

v-T V + 1 
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The second term on the L.H.S. of (6) can thus be written as; 

(7) Pr{C„ <C„; Vi/ € L„) = £ Pr{C m -< C„; Vw € L m \ 

Vm 6 L n 

V-2 i 

+ L Z —^ Pr l C '" = c « < c « ; Vw r e L". v , Vw e i" v ) 

v-T v + 1 

= I V- 

V/w € 

Substituting this into (6) gives 

+ I = I 

v »>«/.,, 

V 

i.e., V i n m = 1; and this proves Theorem 1. 


REMARK: If = 1 for any n € L then n is an absorbing state and the solution to 
Equation (4) is trivial. However, the occurrence of such a situation in real life would be quite 
rare —we therefore assume for our present problem that 0 < t„ m < 1, V(«,m) € L 

LEMMA I. If = 0, («,,/)>) € L. for any n t ^ n 2 and t„ „ ^ 1, then: 

(a) either t n ^„ 2 = 0, V«, € Z., and n 2 is a transient state which can be deleted from the 
chain, or 

(b) there exists at least one /t 4 ^ «, such that i„ > 0 and therefore state n 2 can be 
reached from state rt\ via state n A , hence, all n € L form an irreducible markov chain 
(Note that « 4 € L). 

PROOF. Lemma 1(a) is a compliment of Remark above and does not require a proof. 
For Lemma Kb), if there exists « 4 ^ such that i n ,,, > 0. then there also exists at least 


one A' tuple (/| ,i 2 ./\) such that <T(« 4 ,/ 4 ) > C’Ui 2 .i 2 ) ^ C(« 5 ,/,) ^ 

C'(»|,/|) < C(n k ,i k ): VAr > 4. V/ ( f C, I < j < N. 

Further let 

(vii) L q A {/1 = 0; V« € L ) 

(viii) L n A L - L n 

<ix) L n „ 4 Z, 0 n Z.„ and let N 0 be the number of elements in L 0 . 

Define an 7V 0 x 7V 0 matrix T 0 = (i nm ), V(n,m) 6 Z. 0 < and also define an An vector 
V; - f^f); Equation (4) can now be stated in a modified form, for only the epochs 


"niamed in Z.,,. as; 

‘' A ( f 1 = Ax T(j. 

I MM \ 2 i,, > (I, Vn f To, hence all n € Z» are aperiodic states. 


MHt k wh: 


NAVAL RESEARCH LOGISTICS QUARTERLY 




480 


\ S \1IA 


PROOF: From the cost structure, it is apparent that if r„ „ = 0 then r m „ -= 0, Vm 6 L , 
and hence n_Q T n i.e., state n would have been deleted from the chain considered in equation 

(8) (i.e., in Z 0 ). Therefore, t„„ > 0, Vn € Lo, and hence Lemma 2 follows. 

THEOREM 2: There exists a limiting probability vector A 0 = lim Ao'such that 

d~ °o 

(9) A 0 = Aq x T 0 

has a positive solution which is unique. 

PROOF: The necessary and sufficient conditions for the existence of a positive and 
unique solution A 0 is that the chain should be irreducible and aperiodic. Both conditions were 
established by Lemma 1(b) and Lemma 2 respectively. 

The interest of this paper is in LT s choice of arrival time in the long run, A 0 , which can 
now be solved for in Equation (9). 

2.2.2. A Special Case 

In the general model it was assumed that if U arrived at epoch n on the cf h day, he would 
not change his arrival time to any other epoch if C„ ^ C m , Vm 6 L n . Suppose we modify this 
assumption and now assume that when C„ = C m for any m ^ n we would not rule out the pos¬ 
sibility of V changing his arrival time to epoch m. In that case we shall attach equal chances of 
U changing his arrv;.; time to epoch m and to him not changing his arrival time from epoch n. 
Thus, l/s strategy will now be modified such that the transition probabilities that define his 
strategy will be given by: 

(3a) t nm = Pr{C m < C u ; Vu € L m \ + 

V —^—T PriCm = C m < C H ; Vm, € L m 

v + 1 

Vw 6 L m v }: Vn,m 6 L, 

where 

(x) L m y A 1/1/ = *n|,w 2 , . • • , m v ) such that 

c m = c m| = C„ 2 = ... = C„ v ; 1 < v < A - 1. and 
Z. m , v c L . and 

(xi) v A L m - L m v . 


It is immediately apparent that the right hand side of Equation (3a) is independent of n, 
hence t nm = r v m ; V(n.m,\ ) € L. Let t m A t„ m \ V(n.m) € L. For this special case, therefore, 
all the rows of T are identical, and the steady state solution is given by 


( 10 ) 


a„ 


!„ ; Vn € Lt) 
0 Vn 6 L 0 . 
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2.2.3. Illustrative Example 


For illustrative purpose let, 

/ = 2. N = 3. r = 3. M = 1, 
A„ = 0.3, A| = 0.5, A, = 1.0, 


Let; 


(ID 


C(n.i) 


C u x (/) + 


C] x (r — a — i) . 
C / * (// + /'— T ) 


and let (,. = 1.0. C, = 1.0, = 2.0. 


n + i ^ r 
n + i ^ t 


By using results from Minh [4] we obtain the distribution of queue lengths and from Alfa 
[1], the distribution of waiting is obtained as 


W’ 



_ " _ _ J 


1 

2 


0 

0.73 

0.60 

0.57 

1 

0.24 

0.36 

0.39 

2 

0.03 

0.04 

0.04 


Hence, F" — r *= 1,2,3. 

The total cost C(n,i) are given as 


C(n,i) 



n 


jj 

2 

3 

1 

3 

2 

3 

2 

4 

5 

6 

3 

7 

8 

! _2_ 1 


The transition probabilities are thus given as 

(12) = V\ (Fj + Fj ] + Fj [Fj + Fj] [Fj + F j) + Fj FjFj. 

(13) (| j = Fj + FjFj [Fj -f Fj]. 

(14) r,F j [ F 2 ' + F,'] [ Fj + Fj] + FjFjFj, 

(15) / 2 , = F, 1 [Fj + Fj] fF 2 + Fj] + F 2 ' [F 2 2 + F 3 2 ] [F 2 3 + Fj) 

+ FjFjFj + F,'[F 2 2 + F, 2 1 Fj/2, 

(16) r 2 . 2 = F, 2 + Fj Fj l Fj + Fjl. 

(17) t 2 .., - FjFj [Fj + Fj I + [Fj 4- Fj) Fj Fj + Fj Fj [ Fj + Fj )/2, 
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(18) /, i = [V\ + V} ] (Vj + VlWl + [V} + F, 1 I (F, 2 + FflF/ 

+ V}VjV}. 

(14) , 3 2 = F, 2 + F 2 2 F,' [ F 2 3 + V]\, 

(20) i 7 ,’ [f 2 2 + f, 2 | + f/ j + r,’j.' 2 rJ. 

Thus, this leads to the transition matrix, 

0.333 0.605 0.062 

(21) T = 0.256 0.605 0.139 

0.167 0.605 0.228 

from which we obtain 

(22) A =[0.265, 0.605, 0.130], 

Note that for this example L 0 = L, hence T 0 = T and A 0 = A. 

3. CONCLUSION 

This model, based on its assumptions, can only be used to predict the customer’s arrival 
time at the system when all other customers previously using the system have settled for a par¬ 
ticular arrival time and will not change it under any circumstances. However, when all the 
other customers change their arrival times as a result of the new customer's "interference" then 
there will be a slight modification to the problem. 

If other customers change their arrival times everyday, then the distribution of W" will 
change from day to day and, hence, T and T 0 will also change from day to day. If we let T d or 
represent the transition matrix for the strategy for the (d + l) th day then If s choice of 
arrival time can be reformulated as 

(23) A rf+I = \ d x T d . 

However, with (23) existence of a steady state solution cannot be assured. If we further let 
there be more than one server in parallel such that a customer develops a strategy noi just for 
selecting his arrival time only but for joint selection of arrival time and of the server to serve 
him, then the problem becomes similar to that in Alfa [2], 
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ABSTRACT 

This paper develops a forward algorithm and planning horizon procedures 
for an important machine replacement model where it is assumed that the tech¬ 
nological environment is improving over time and that the machine-in-use can 
be replaced by any of the several different kinds of machines available at that 
time The set of replacement alternatives may include (i) new machines with 
different types of technologies such as labor- and capital- intensive, (ii) used 
machines. <iii> repairs and/or improvements which affect the performance 
characteristics of the existing machine, and so forth 

The forward dynamic programming algorithm in the paper can be used to 
solve a finite horizon problem The planning horizon results give a procedure 
to identify the forecast horizon T such that the optimal replacement decision 
for the first machine based on the forecast of machine technology until period 
T remains optimal for any problem with horizon longer than T and. for that 
matter, for the infinite horizon problem A flow chart and a numerical exam¬ 
ple have been included to illustrate the algorithm 


1. INTRODUCTION 

In our previous paper [5l, we developed forward algorithm and planning horizon pro¬ 
cedures for a machine replacement model under the assumption that only one kind of machine 
is available for replacement in any given period. We showed that there exists a forecast horizon 
T, such that the optimal replacement decision for the first period (based on the forecast of 
machine technology until period T remains optimal for any longer (than T) horizon and, for 
that matter, the infinite horizon. 

In this paper, we relax the assumption of a single possible replacement alternative by mul¬ 
tiple alternatives. That is, we develop a model in which the machine-in-use can be replaced by 
any of the several different kinds of machines available at that time This model can deal with 


VOL. 29, NO. 3, SEPTEMBER 1982 


483 


NAVAL RESEARCH LOGISTICS QUARTERLY 





484 


sell \NI> \\l) S SI nil 


many realistic situations. For example, the set or replacement alternatives may include (i> new 
machines with different types of technologies such as labor- and capital- intensive, (ii) used 
machines, (iii) repairs and/or improvements which affect the performance characteristics of the 
existing machine, and so forth. It should be noted that the set of replacement alternatives can 
change from period to period. 

This multiple alternative replacement model represents a major generalization of the vari¬ 
ous single alternative replacement models developed by Terborgh [71, Thompson [81. Sethi and 
Morton [61, Gordon (1). For a survey of the extensive literature on the machine replacement 
models, the reader is directed to Rapp [41; see also Pierskalla and Voelker [31. 

In the next section, we formulate the multiple alternative replacement model under 
improving technologies. In section 3, we develop an efficient forward algorithm for solving 
finite horizon problems. Section 4 develops a variant of the regeneration- monotonicity pro¬ 
perty [2, 5] and planning horizon procedures for the model. A flow chart summarizes the com¬ 
plete solution procedure. A numerical example is solved in Section 5 to illustrate the steps of 
the forward algorithm and the application of the planning horizon procedure. We weaken the 
assumption of improving technologies in Section 6 and show that planning horizon procedures 
can be adapted to this case. Section 7 concludes the paper with some important remarks. 

2. MODEL FORMULATION 

Consider the situation of a production shop which must keep a single machine of a partic¬ 
ular capacity at all times To run this machine, the shop incurs operating expenses which may 
include labor cost, electricity cost, maintenance cost, depreciation, and so forth. Usually, the 
performance of the machine deteriorates, i.e., operating expenses increase over time, so that 
the shop might consider selling the existing machine for its salvage value and buying a new 
one. This new machine is selected from various different alternative machines available If a 
major repair is decided on the existing machine, then this alternative can be considered as if a 
new machine of another technology is purchased at the cost of repair plus the salvage value of 
the existing machine. Once the new machine is purchased, the same situation repeats with the 
new machine. Consequently, for any finite or infinite horizon, the shop will make a chain of 
replacement decisions [4], (51. The problem of the shop, of course, is to find simultaneously 
the optimal times of these replacements and the types of machines selected at these times. 
Obviously, these decisions will depend on the prices of the future machines and their per 
period operating expenses over time. More precisely, we are considering the following model 
which the shop must solve: 

Let A * denote the machine of technology h available at time t. Let there be machines of 
N alternative technologies available to choose from in any period. It is. noted that the extension 
to the case, when the set of alternative technologies is changing over time, is straightforward 
and will be described in Section 7. 

Let 0, h k , k ^ r, be the operating cost of the machine A, h in period k. With reference to 
our previous paper [51, we note that 

O* = tt/ 1 - S, h , + M?, 
and 

0,\ - - S* k + 
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where it/' is ihe purchase price (il is assumed thai the new machine A, 1 ' is purchased at the 
beginning of period /) and S,'' k is the salvage value at the end of period A of the machine A". 
T he term A//' t denotes all the costs, which are specific to the technology h, of operating the 
machine A 1 ,' in period A excluding the cost of depreciation (or loss in salvage value from period 
A - I to A). 

To obtain a forward algorithm, we must solve the multiple alternative replacement model 
for any given finite horizon P. Suppose the shop uses n machines during the interval <1 . P>. 

Let t 1 . h . i„- t„ be their salvage times and let h\, h 2 . h„ be their technologies. 

Note that r, ^ l for i = 1, 2, ... , (n — 1) and t„ = T, since we assume that the shop goes out 
of business at the end of the given finite horizon P. Define to = 0. Also note that 

<p,q> = [p, p + I. q] with nonnegative integers p and q with q > p. We can now state 

the finite horizon problem as 

n - 1 '/.l . 

(U Minimize £ £ Of'+Y.k 

,-n + \ 

«. »t. h . >n -1 

h], h 2 . ■ ■ ■ , h„. 

In the next section, we develop a forward dynamic programming algorithm to solve (1). 

3. FORWARD ALGORITHM 


We start with defining the following terms: 

Purchase Point. A period t is defined to be a purchase point (or P-poinii if a machine is 
purchased at the beginning of period t. 

Regeneration Point A period t is a regeneration point (or R-point) if a machine is sal¬ 
vaged at the end of period t. 

It is easy to see that any P-point is immediately preceded by an P-poinl. This property 
allows us to develop an efficient forward algorithm [2]. 


We need to introduce the following notation: 


cm 

CHj.T) 


C*( T) 


(3) 

Cj(T) 


— minimum cost of the T-period problem (1), 

= the totai operating cost for the machine 
/f/Vi in the interval <j + I. T> 

= £ o; + ,* 

k-j* 1 

= minimum total cost of the P-period problem when (./ + 1) is the last 
P-point and h is the technology of the last machine 

= C(J) + C h (j T) 

= minimum total cost of the P-period problem with <J + 1) as the last 
purchase period, 
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Clearly, 

(4) C,m = min CHT) 

h 

and 


(5) C(T)= min CAT) 


Equations (2)-(5) complete the statement of the forward algorithm. It is possible to sim¬ 
plify the computation of C*(D; using (2) in (3) gives the following sequential procedure 


( 6 ) 


C?< T) = 


C*(F — 1) + o; + ,. r for j < T - 1 

c(r- i) + of r j = r- r 


4. PLANNING HORIZON RESULTS 

To obtain planning horizon results, we need to assume that every technology is an 
improving technology. Thus, for technology A, we assume that 

(7) 0/L u > 0** for k Ss ./ + I. 

In other words, we assume that the operating cost in a period for a machine of technology A is 
lower than the operating cost for an older machine of the same technology, except perhaps in 
the first period of operation of the new machine when the amount of depreciation on the new 
machine is usually high. 

Lower Bound Monotonicity Property: 

Let j*(T) denote the latest next to the last /{-point in the optimal solution of the T- 
period subproblem; note that period T is the last /{-point. By the tower bound monotonicity pro¬ 
perty we mean that the lower bound of j*(T) increases monotonically with T. To prove this 
property, we must also define j* h (T ) to be the latest next to the last /{-point for the optimal 
solution of the C h ( 7')-problem, where the C h i D-problem is the T- period problem subject to 
the constraint that the last machine be machine of technology h. 

THEOREM 1: Under the assumption (7) of improving technology over time, j* h (T) 
satisfies the regeneration monotonicity property, i.e., 

(8) j* h (T + 1) ^ j'HT). 

PROOF; Let j* h (T + I) » a and j* h (T) — b. It is easy to see that C^iT) > C£(T). 
Furthermore, if a < A then 0* +) r+ , ^ 0*+ 1 . r +1 • 

If a < A, then using (6) we can write 

C*(f+ 1) - C*( T) + 0 0 \ i r+l 
^ CUT) + 0»\,. r+1 
- CZ(T + 1). 
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This contradicts the optimality of the C^(T + 11-solution, implying that a ^ b. This completes 
the proof. 

THEOREM 2: If we define 

(9) X(D = min (/**(D), 

h 

then A(D represents a lower bound of j*( T 4- 1) and X(D is monotonically increasing in T, 
i.e., 

(10) j *( T + 1) Js X(D Si X(T - 1). 

PROOF: The proof easily follows from (8) and (9). 

Regeneration Set: 

The regeneration set of the C(T )-problem can now be defined as 

on r(t) - set - {xm, xm + i.r-il. 

We note that a set is known as an R(T )-set if it contains at least one /{-point of an optimal 
solution to any problem with horizon T or longer. 

In many cases, it is possible to reduce the size of the /?(D-set. This reduction makes it 
easier to obtain planning horizons as can be seen in the planning horizon theorem below. The 
reduction procedure requires computations of the S'TD-sets for all h, where an S h (T )-set is 
defined to be a set of regeneration points such that j* h ( T + /) € 5'’( D-set for all periods 
(T + I) with j* h (T + /) < T and / > 0. It is convenient to express this procedure by the flow 
chart on the next page. 

The reduced R ( D-set can be obtained as follows: 

(12) R (T )-set = U S h ( D-set. 

h 

It should be noted that the R( D-set in (12) cannot be bigger than the R( D-set in (11). We 
now prove the following important result for the set 5 found in the flow chart. 

THEOREM 3: If a period t £ the 5-set found in the flow chart for technology It. then 
./**( T + /) ^ t for any value of / > 0. 

PROOF: If t 9 the 5-set found in the above flow chart, then there is a period it. 
T - I > u > t, such that C*( D < C,*(D. From the improving technology assumption, we 
also have 0//+ u < 0,+ j * for k > T + 1. We can now write for any period T + /: 

T+l 

c„*(r + /) - cm T) + £ 0* +l .* 

k-T +1 
T+l 

< C,*(D + I 0,* +1 ,* 

*-r+i 

- C, h (T+ I), 

implying that j mh (T + I) r* t. This completes the proof. 
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Flow Chari for Computing the S*< D-sel. 


Planning Horizon Theorem: 

Let f*(r) denote the period when the first machine is salvaged in an optimal (re¬ 
solution. We now state the planning horizon theorem which gives us a stopping rule for the 
forward algorithm to find the optimal replacement time in an infinite horizon problem. 

THEOREM 4 (Planning Horizon Theorem): If ,r(r) = t* = constant for all r € R(T)- 
set (12) for some T such that j*(T) ^ 0. then the salvage lime of the first machine in an 
optimal infinite horizon solution is r*. Furthermore, the optimal technology It* is determined 
by the condition 

(13) C,* 1 Of) = C () 0f>. 


The first replacement time /* is called the planning horizon and T is called the forecast 
horizon. We need only to forecast the relevant technology up to period T to find the optimal 
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technology h\ and the optimal replacement time /* for the first machine. The reader is referred 
to [2] for a discussion of such theorems. 

We now present a flow chart for the complete forward algorithm along with the planning 
horizon procedure. We will use this flow chart to solve a numerical example in the next sec¬ 
tion. 



Flow Chart for Implementing the Planning Horizon Procedure 


VOL. 29, NO. 3, SEPTEMBER 1982 


NAVAL RESEARCH LOGISTICS QUARTERLY 














490 


S ( HAND AM) S SI III! 


5. NUMERICAL EXAMPLE 


We now illustrate the planning horizon procedure by solving a numerical example. This 
example assumes that machines of two different technologies are available for replacement at 
any time. The operating costs O ,J and O, 2 , are as shown below: 

O,] 0, : , 


j 


160 120 

130 

140 

150 

160 

160 

120 

130 

140 

150 


140 

105 

115 

125 



140 

105 

115 




120 

80 





120 


j 


300 60 

70 

80 

90 

100 

300 

60 

70 

80 

90 


200 

50 

60 

70 



400 

30 

40 




400 

30 





200 


Note that technology I can be interpreted as a labor intensive technology and technology 
2 can be interpreted as a capital intensive technology. The calculations of C,HT), C 2 (T), and 
C,(T) are shown below: 



Sample calculations for T = 3 are shown below: 

Cd (3) - Of, + Of 2 + 0,'j = 160 + 120 + 130 - 410, 
C,'(3) = C(l) + 0 2 ' 2 + 0 2 ' 3 = 160 4- 160 + 120 = 440, 
C\ (3) = C(2) + Oh - 280 + 140 = 420, 

7*'(3) = 0. 
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C'S (3) = Or, + Of: + Of, = 300 + 60 + 70 = 430. 

Of (3) = CM I) + Of, + 0 2 2 , =• 160 + 300 + 60 - 520, 

CM 2 <3> - CM2) + Of, - 280 + 200 = 480, 

/*-(3) = 0, x(3) = min{y* 1 (3). ,/* 2 (3)} = 0, 

C„(3) = mint0,1 (3). C’,r (3)) = 410. 

0,(3) = minIC', 1 (3), 0, 2 (3)) = 440, 

C M3) = mini CM 1 (3). <T? (3)) = 420. 

7*13) = 0, /* (3) = 3. 

The planning horizon theorem is not satisfied for T ^ 5. For example, for T = 5. we 
have S'(5) = (4), S 2 (5) = {2,3,4}, and R( 5) = {2,3,4}. Since f*(2) ^ /*(3) ^ /*(5). the 
planning horizon theorem is not satisfied. For T = 6, we have s'(6) = (5), S 2 C6) = {2.5}. and 

R( 5) = (2,5). /*(2) = /*(5) = 2, so the planning horizon is 2 periods and the forecast hor¬ 

izon is 6 periods. 

6. RELAXING THE IMPROVING TECHNOLOGY ASSUMPTION 

Let /(/?) ^ 1 be a number associated with technology h such that 
0,'L|. J+ * 0* J+k for k > l(h). 

This condition is likely to be satisfied by most improving technologies. The following 
result holds for this assumption: 

THEOREM 5: Let i* h (T) < T — r(/t) (0 if T ^ r(/;)) be the latest period such that 
C] */;,,, < C,(T) for all j € <0, T — i(h)>. then ./**( T + 1) > i* h (T). 

PROOF: Assume to the contrary that j* h (T + I ) — a < i' h (T ) It is easy to see that 
O a \, r .i 5s O," V\(T) + 1. T + 1 and C h a (T) > C* */t (ri m. We can now write 

C */t | t ,(T + 1) - C* */t, r ,(r) + 0* • A (n+ , r+ i 

=$ c a ft m + o* +l . r ,, 

= c fl "(r + n. 

implying that j* h (T + 1) ^ i* h (T). 

Regeneration Set. 

Let k(T) = min [/**(7~> 1, then the regeneration set of the CMD-problem can be defined 

h 

as 

R(T) — set HUfUm + l. T - 1}. 

As in Section 4, it is possible to find a reduced regeneration set. For this we first find the 
S h ( T )-set as below. The reduced regeneration set can be found by using (13). 
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Flow Chart lor Computing the S"( 7T-sot 


7. EXTENSIONS AND CONCLUDING REMARKS 

In this paper we have solved the infinite horizon replacement model with multiple alterna¬ 
tives. We have obtained planning and forecast horizons for this model under reasonably gen¬ 
eral conditions Furthermore, the model can be easily adapted to situations where the set of 
alternative technologies is changing over time. This is done by setting Of], = for the ima- 
i ginary A]' machine for all / < t for a technology which appears at lime i for the first time. If a 

* technology h disappears in period t. then we let O'', = <» tor fictitious machines A]', for all 

,/' ^ i With these definitions, the model developed in the paper is applicable to the situation of 
changing sets of technologies. 

Another situation to which the model can be easily adapted is the situation in which a 

switch over cost k ah occurs whenever a machine of technology a is replaced with a machine of 

technology b. For this, we do require the assumption that k ab is separable, i.e.. 

k uh = k“ + k h . 

In this case, we can adjust the salvage value of the existing machine downward b\ amount k" 
and adjust the price of the new machine upward b\ amount k h . 

Finally, we must slate that the extensions of the machine rcpi.ucmcni model m ..m r. . 
ous paper [ 51 can also be solved in the multiple alternative case m a t.uiU stigirt 'w 
manner 
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FLOWSHOP/NO-IDLE OR NO-WAIT SCHEDULING 
TO MINIMIZE THE SUM OF COMPLETION TIMES 


I. Adiri and D. Pohoryles 
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Haifa, Israel 

ABSTRACT 

This paper deals with flowshop/sum of completion times scheduling prob¬ 
lems, working under a "no-idle" or a "no-wait" constraint, the former prescribes 
for the machines to work continuously without idle intervals and the latter for 
the jobs to be processed continuously without waiting times between consecu¬ 
tive machines. Under either of the constraints the problem is unary NP- 
Complete for two machines. We prove some properties of the optimal 
schedule for n/2/F, no-idle/IC,. For n/m/P, no-idle/X C, and n/m/P, no- 
wait/EC, with an increasing or decreasing series of dominating machines, we 
prove theorems that are the basis for polynomial bounded algorithms. All 
theorems are demonstrated numerically. 


INTRODUCTION 

The nonpreemptive, flowshop, sum of completion-times scheduling problem is: n jobs 

U\,Ji . J„) have to be processed by m machines (M U M 2 , .... M m ). Job 7 j( 

i “ 1,2.n, consists of, at most, m operations (0,i,0, 2 , ... , 0 fm ). Operation 0 u which 

precedes Q (J+1 , has to be processed uninterrupted for time units, on M r j - 1,2. m. t 0 

is a nonnegative integer — r /y * 0 if 0 tj is missing and positive if it exists.* Two operations of 
the same job cannot be processed simultaneously and a machine may process at most one job at 
a time. Find the operation sequence on each machine, that obeys the problem constraints and 
minimizes the sum of completeion times. The problem is designated n/m/F/EC*,+ where F 
stands for flowshop discipline and C t is the completion time of job J,. 

The no-idle constraint: machine M k , k - 1,2, ..., m, works continuously without idle 
intervals. 

The no-wait constraint: job J h i - 1,2, ... , n is processed continuously without waiting times 
between consecutive machines. 

Both constraints arise in real life situations. Examples of such scheduling are: (i) under a 
no-idle constraint — use of very expensive equipment (.e.g., a computer and its peripheral dev¬ 
ices) with the fee determined by the actual time consumption; (ii) under a no-wait constraint 
— in metal-processing industries (e.g., hot rolling) where delays between operations interfere 
with the technological process (e.g., cooling in the above case). 

"For further elaboration of the meanings and influence of zero processing times see Hefetz and Adiri (41. 
tFor notation and classification of scheduling problems we follow Lenstra [5] and Rinnooy-Kan [61. 
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I AOIRI and d. pohoryles 


While ihe no-wait case has been widely discussed in the professional literature, fl) (61, 
the no-idle case is defined here for the first time. This paper deals with flowshop disciplines 
with a no-idle or a no-wait constraint where the objective function being minimum sum of 
completion times. 

Garey, Johnson and Sethi [3] proved that »/2/F/XC, is unary NP-Complete by construct¬ 
ing an instance that is no-idle and no-wait, and thus proved at the same time that of n/2/F, 
no-idle/XC, and n/2/F, no-wait/LC,. Section 2 deals with some properties of the optimal 
schedule for n/2/F, no-idle/XC,. Sections 3 and 4 are devoted to proofs of theorems that 
underlie polynomial bounded algorithms for n/m/P , no-idle/XC, and n/m/P, no-wait/XC, with 
an increasing or decreasing series of dominating machines, respectively. 

1. COMPLEXITY OF FLOWSHOP/NO-IDLE OR NO-WAIT/SUM OF 

COMPLETION TIMES PROBLEMS 

Garey, Johnson and Sethi (3] proved that 3-partition is reducible to n/2/F/T.C,. How¬ 
ever, as the constructed instance of n/2/F/XC, happens to be no-idle and no-wait, thus at the 
same time proved the unary NP-Completeness of n/2/F/ZC h n/2/F, no-idle/XC, and n/2/F, 
no-wait/EC,. Moreover, only minor modification of the proof is needed for proving the NP- 
Completeness of the first two problems where missing operations are prohibited. Specifically, 
replacement of zero processing times on M\ (missing operations) by c > 0 (existing operations 
with infinitesimal processing times) and shifting of all c to the beginning of the schedule on 
A/], proves the unary NP-Completeness of n/2/F, i tj > 0/XC, and n/2/F, t u > 0, no-idle/XC,. 
The complexity of n/m/F, t 0 > 0, no-wait/XC* for fixed m ^ 2 is an open problem. 

2. CONSTRUCTION OF A NO-IDLE SCHEDULE 

We distinguish two ways of constructing a no-idie schedule. Both are implemented con¬ 
secutively on the machines, starting with the second, proceeding to the third and so on until 
the last machine. 

(i) Right shifts. We shift to the right every operation that precedes an idle interval 
(maintaining the constraints of the problem on M k , k > 2, further right shifts might be 
needed), until all idle intervals have been eliminated. 

(ii) Left shifts. At first we schedule the operations on M k , k > 2 in a single stretch (a 
block) without idle intervals, starting when the last job on A/*_i has been completed. After¬ 
wards we apply a maximum left shift (without violating the constraints) to the whole block. 

Let us define a blocking job on Af k , k > 2, in a no-idle schedule as the first job on M k 
that prohibits shifting to the left. 

3. PROPERTIES OF THE OPTIMAL SCHEDULE FOR n/2/F, NO-IDLE/XC, 

Conway, Maxwell and Miller [2] fundamental theorem for flowshop scheduling that states 
that an optimal schedule exists for a n/m/F/B problem (8—any regular measure of perfor¬ 
mance) with the same processing order on the first two machines [2, p. 81], holds for the case 
under discussion. Thus, the set of permutation schedules for n/2/F, no-idle/XC, is a dominant 
one. 
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THEOREM 1: If in the optimal schedule for n/2/F, no-idle/EC,, the blocking job, is the 
last one, then the optimal schedule is according to SPT (Shortest Processing Time) on M 2 , 
except for the blocking job that is the one with the minimum processing time on M 2 . 

PROOF: Let us denote, a, - / (l , 6, - t l2 , S ilt starting time of J, on M k \ and the square 
brackets indicating the place of a job in the sequence, for example, Sm 2 —starting time of the 
first job on M 2 . Since the blocking job is the last one, we have (see Figure 1). 

W %2 " X a l<l ~ L *I<1 “ K + *(»l" 

/-i i-T 

where 

n n 

K - £ a, — £ b, — Constant. 

“ i-\ 



Figure 1. Two machines no-idle schedule where the 
blocking job is the last. 


The sum of completion times takes the form 

( 2 ) £ Ci - S[|] 2 + Am + 5 (i] 2 + *lil + *121 + • • • + 5(1)2 + *[/) 

- «5(1|2 + (» - i + 1 )*[/)■ 

Substitution of (1) yields 

(3) £ C) “ nK + nb( B | + (n — / + 1)A[/| 

it -1 

— nK + (n + 1)*1b1 + ^ (n — / + 1)*(#)- 

To minimize (3) />(„] should be the smallest and &(,), / — 1.2. n — 1, a nondecreasing 

sequence, thus the optimal schedule is given by the sequence 

*1*1 ^ *11) < *121 < - • • < *tir—11- LI 

THEOREM 2: The jobs that (i) precede (ii) succeed the blocking job in the optimal 
schedule for n/2/F, no-idle/1 C, are ordered according to SPT on M 2 , provided the no-idle con¬ 
straint is not violated. 

PROOF: Let us assume that the blocking job in the optimal schedule is the «r-th in the 
sequence. 
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I ADIRI AND I). POIIORYLES 


Mi 




M? 



Figure 2. Optimal schedule for n/2/F, no-idle/EC, where the blocking job is 3|,|. 


w 5 ui2 “ x a i«i ~ *i>i 

The situation here is rather simitar to that of Theorem l with w replacing n. Thus, in the 
optimal schedule we have Am < A( 2 | < ... < *[»-n- 

(ii) For / - n + 1, ir + 2, .... n, we have 

C [<1 " C lirl + X *(/!• 

j~rr +1 

Thus, 

X Cli) ” ( n ~ w )Qirl + X ~ 

/“ 1T + 1 /—W+l 

We have that in the optimal schedule the jobs that succeed the blocking job obey the SPT rule 
on M 2 , provided the no-idle constraint is not violated, A[ w+ n < A[„+ 2 l < .. • < £(*]. □ 

4. n/m/P, NO-IDLE/I C„ m> 2, WITH A SERIES OF DOMINATING MACHINES 

Machine M k dominates M, if 

(4) min t ik > max 

/ / 

In abbreviated notation 

M k > M r . 

A scheduling problem with a series of increasing [decreasing] dominating machines is one 
where M\ < M 2 <... < M m lM\ > M 2 > ...> M„\. 

For n/m/F\ no-idle/IC ( , m > 2, the set of all permutation schedules is not a dominate 
set. However, for the sake of solvability (development of polynomial bounded algorithms for 
special cases) we confine our search for optimality to permutation schedules, and the problem is 
designated n/m/P, no-idle/EC/. 

Note that for m - 2 the set of permutation schedules is a dominate set. 

Let y, k be the processing time of the Mh job on M k where the order is according to SPT 
on the last machine, M m . 


NAVAL RESEARCH LOGISTICS QUARTERLY 


VOL. 29, NO. 3, SEPTEMBER 1982 



SCHEDULING TO MINIMIZE SUM OF COMPLETION TIMES 


499 


THEOREM 3: The optimal permutation schedule for n/m/P, no-idle/IC, with an 
increasing series of dominating machines* is according to SPT on the last machine, except for 
the first job that satisfies 

1 m—1 

v i -" z 

*-l 


j -1 


yjk + O' - Dyym - ZJ, y»i 


PROOF: For this case we have 


(5) 


c, - «5(i) m + (« - / + 


where S[i] m is the starting time of the first job on the last machine, and 
(6) 5[i| m - /[,],. 

Thus, for a given first job the optimal permutation schedule is SPT on M„, 

'I2)m < '(31 m < ••• < 

Selecting the >th job in the SPT sequence on M m to be the first, we have 

O) Z C ‘ ” " Z y Ji + 0 - - Z y>« + Z 

Since the last term in (7) is not affected by the choice of the first job, the latter is taken so as 
to satisfy 


m -1 j—\ 

V i “ « Z * + 0 “ 1) y Jm ~ Z T/m 

*-l f-l 


□ 


EXAMPLE 1: A 5/A/P, no-idle/IC, with an increasing series of dominating machines 
with processing times as per Table 1. 


TABLE 1 — Processing Times 


Jobs 

Machines'"^ 

J\ 

Ji 



h 

Mi 

~T~ 

_ 

~r~ 

T 

~7 

Mi 

4 

5 

6 

4 

6 

Mi 

7 

9 

8 

8 

9 

a/ 4 

ll 

13 

14 

10 

12 


A / 4 > A /3 > M 2 > Mi, thus we have an increasing series of dominating machines. 
SPT sequence on A / 4 is 4-1-5-2-3. 


'Since M k+I > M k , k - 1,2.m - 1, the no-idle constraint is not effective and for the case under discussion, the 

two problems n/m/P, no-idle/IC, and n/m/P/ZC , are equivalent. 
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The firs! job, say J„ satisfies P, — Min Vj, is the Ath in the SPT sequence on A/ 4 . 


V, - 5(2+44-8) . - 70 

V 2 - 5(3+4 + 7) + 11-10 - 71 

y 3 - 5(2+6+9) + 212 — (11 + 12) -88 

V 4 - 5(1 + 5+9) + 3-13 - (10+11 + 12) -81 


y s - 5(3 + 8+8) + 4 14 — (10+11 + 12 + 13) - 95 


The job chosen to be the first is the first in SPT sequence on M 4 (F, - min F y ), namely J 4 . 

j 

The other jobs are scheduled according to SPT on M 4 Thus, the optimal permutations 
schedule is 4-1-5-2-3 with T C, — 240, and takes the form as per Figure 3. 



THEOREM 4: The optimal permutation schedule for n/m/P, no-idle/IC, with a decreas¬ 
ing series of dominating machines is according to SPT on the last machine except for the last 
job that satisifies 

mini Vj - n T y Jk - (n - j)y, m + T y im 
J l *-2 i-y+i 

PROOF: Since M\ > M 2 > ... > M m the blocking job on M k , k — 2,3.m, is the last 

one, we have 

m—1 

(8) S[ii m - (5(i](*+n ~ ■S’tiifc) 

«l-l [ PI H-l 

"II ^1/1* ~ X *I<K*+I) 

*-I (/-T /-I 

“ X ~ 7*+1 + J i<iiu+i)) — K + 2* 

*-l k-l 


where T k is the total processing time demanded on M k , T k - constant 
K-T,~ T m . 


n 
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Substitution of (8) in (5) yields 

<9) £ C, - nK + n £ / (b |(* + d +£(«-/+ 

i-l L-l /-I 


Given the job in the last (n-th) place, we have 

(10) jj c, = Af, + jj (!»-/• + i)/ w „ 

m—l 

where A - ! — constant - nK + n £ /(„]<*+,, + 

Minimum of (10) is obtained by an SPT sequence on M m , 

'Him < t\2\m < • • • < 

Choosing the >th job in the SPT sequence on M m for the last place, we have 

(H) T C, - nK + n £ y y <*+i> - (n - + T (n - i + l)y /m . 

/—T *-i i-j +1 i-i 

The last term on the right hand side of (11) is independent of the choice for the last place, thus 
the job is taken so as to minimize 

| y j-» X yjoc+n - jhjm + D 


COROLLARY 1: The optimal schedule for n/2/F , no-idle/IC, where M\ > M 2 is 
according to SPT on M 2 except for the last job that is the one with minimum processing time 
on Mi- 


PROOF: As was pointed out the two problems n/2/P , no-idle/IC, and n/2/F, no- 
idle/IC,, are equivalent. Substitution of m » 2 in Theorem 4 yields 

v \ “ min I Vj - jy n + T yd 
1 \ '-/+' ) 

Thus, the last job is the first in SPT sequence on M 2 , This result is in agreement with 
Theorem 1—since A/j > M 2 the blocking job is the last one and should be the smallest. □ 

EXAMPLE 2: A 5/4 /P, no-idle/IC; with a decreasing series of dominating machines 
with processing times as per Table 2. 


TABLE 2 — Processing Times 


"N. Jobs 

Machine^ 

A 

A 

A 

A 

A 

M\ 

Ti" 

TF 

' 14 

To' 

Tr 

M 2 

7 

9 

8 

8 

9 

Mi 

4 

5 

6 

4 

6 

M t 

3 

1 

3 

2 

2 
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A/i > M 2 > M } > M 4 , thus we have a decreasing series of dominating machines. 
SPT sequence on A/ 4 is 2-4-5-1-3. 

The last job, say J h satisfies V t = Min Vj, is the Ath in the SPT sequence on M 4 . 

j 

V y = 5(9 + 5 + l) - 4-1 + 2 + 2 + 3 + 3 = 81 


V 2 = 5(8+44-2) - 3-2 + 2 + 3 + 3 =72 

V] = 5(94-64-2) — 2*2 + 34-3 = 83 

y 4 = 5(7+4 + 3) -1-3 + 3 =70 

K 5 = 5(8+6 + 3) = 85 


The job chosen to be the last is the fourth in SPT sequence on M 4 , ( V 4 — min V,), namely J\. 

j 

The other jobs are scheduled according to SPT on M 4 . 


Thus, the optimal permutation schedule is 2-4-5-3-1 with 
per Figure 4. 


s 


V' 


343 and takes the form as 



5. n/ m/P, NO-WAIT/EC,, m > 2, WITH A SERIES OF DOMINATING MACHINES 

We recall that "no-idle" constraint prescribes for the machines to work continuously 
without idle intervals, while the no-wait constraint prescribes for the jobs to be processed con¬ 
tinuously without waiting times between consecutive machines. 

The set of all permutation schedules for an "F, no-wait" problem (t tJ > 0, missing operations 
are allowed) is not a dominant set. However, we confine our search for optimality to permuta¬ 
tion schedules, namely, the problems under discussion in this section, are particular cases of 
n/m/P, no-wait/EC/. Note that for t tJ > 0 (missing operations are not allowed) a feasible solu¬ 
tion for an "F, no-wait" problem is a permutation one, thus "F, ty > 0, no-wait" and "P, t 0 > 0, 
no-wait" are equivalent. As was previously pointed out, n/l/F, no-wait/EC, is NP-Complete 
while the compexity of n/m/F, no-wait, f y > 0/EC,, m > 2, is an open problem. 
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THEOREM 5: For an increasing series of dominating machines the two problems n/m/P, 
no-wait/EC, and n/m/P, no-idle/EC/ have the same optimal sequence. 

PROOF: For n/m/P, no-wait/EC, with an increasing series of dominating machines, we 

have 

C (l | - £ / (1I/ + £ 

/-I j -1 

It b reddily shown that (12) leads to (5) and (6). □ 


A direct consequence of Theorems 5 and 3 is that the optimal schedule for n/m/P, no¬ 
wait/EC, with an increasing series of dominating machines is SPT on M m except for the first 
job that satisfies 

[ m- 1 j-\ 

Mini Vj - n £ y jk + (J - l)y Jm - £ y im . 

1 I *-i /-I 


EXAMPLE 3: A 5/4/P, no-wait/EC/ with an increasing series of dominating machines 
with processing times as per Table 1 (Example 1). 



The calculations and the optimal sequence are the same as in Example 1, but the resulting 
schedule is as per Figure 5 (compare with Figure 3). 

THEOREM 6: The optimal schedule for n/m/P, no-wait/EC,, m > 2, with a decreasing 
series of dominating machines* is according to SPT on M\. 


'Since M k > M k+) , * - 1,2.m - 1, the no-wait constraint is not effective ; for the case si discussion the 

two problems n/m/P, no-wait/IC, and n/m/P/XC, are equivalent. 
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PROOF: Since M k > M k +\, k - l, 2. m - 1, we have for the optimal schedule 

(13) Z C/ “ Z 'in; + 'mi + Z + • • • + Z 'im + Z 'w; 

i-i “ ;-i i-i ;-i 

“ Z Z 'Id; + Z ~ '*'ldl- 

“ j-\ i-i 

The first term on the right-hand side of (13) is a constant, and the sequence 

'mi ^ 'uii < ••• ^ 'in—in minimizes the second term. Morever, r In n > otherwise if 

/(„11 ^ max r(,n the value of the second term can be reduced by interchanging the last job with 

■^l/li* 'OH “ max 't,;i- □ 

EXAMPLE 4: A 5/4/ P, no-wait/IC, with a decreasing series of dominating machines 
with processing times as per Table 2. Since M t > M 2 > My > M 4 we have a decreasing series 
of dominating machines and the optimal sequence is according lo SPT on 4-J-5-2-3, wiih 

£ C, =» 247. The optimal schedule takes the form as per Figure 6. 
i-1 


Mi 

4 | 1 1 5 1 2 I 3 

i 



M~ 

ni m_i s i i 2 i_i 

l-U 

mt 

m 3 

■RKIS 

, LlJ 

m 4 


\ — 3 


24 35 50 

Figure 6. Optimal schedule for Example 4. 
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ON THREE BASIC METHODS FOR SOLVING 
BOTTLENECK TRANSPORTATION PROBLEMS 
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Bonn, West Germany 

ABSTRACT 

For solving transportation problems essentially three types of methods are 
known: primal methods, the Hungarian method and the shortest augmenting 
path method In this paper we present the specialization of these approaches to 
the bottleneck transportation problem and report some computational experi¬ 
ence. 


1. INTRODUCTION 

The classical transportation problem (TP) can be formulated in the following way. There 
are m supply points and n demand points with supply point i capable of supplying amount 

a t U - 1. m) and demand point j having demand bj(J =1, .... n) with £ a, - £ bj. 

< j 

Find the least cost transportation pattern from the supply points to the demand points when the 
cost of transporting a unit amount from / to j is c,j , i.e., 


} 

► 

t 



min H f j/' x u 
<- 1 j- i 

subject to 


(1.1) 

i,*u m 

7-1 

(/ — 1. m ) 


(1.2) 

£ x ij “ bj 

1-1 

0-1. n) 


(1.3) 

\v 

o 

(/ - 1. m\ j = 1, .. 

. , n 


A related problem is the bottleneck transportation problem (BTP). Here a transportation 
time t tJ is specified between each supply point i and each demand point j. Now it is required to 
find a transportation pattern which minimizes the total time necessary for transporting the 
goods from the supply to the demand points. Thus, 

(1.4) min max l/,yUy > 0} subject to (1.1), (1.2) and (1.3) 

Problem (1.4) is often called "time-transportation problem." It occurs in connection with tran¬ 
sportation of perishable goods, with the delivery of emergency supplies or when military units 
are to be sent from their bases to the front. 
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For solving TP we know essentially three types of methods 

— primal methods 

— the "Hungarian” method 

— the shortest augmenting path method. 

Nearly every approach for solving TP can be modified to treat BTP. 

So for instance the primal method of Hammer [10], [10al, for solving BTP is strongly 
related to the method of Klein [11] for solving TP. Starting from a feasible solution both 
methods are looking for "negative cycles" to improve the actual solution. 

In this paper we present these three basic approaches for solving BTP and we report some 
computational experience with these methods. 

2. THE HUNGARIAN METHOD 

This method for solving BTP was first proposed by Garfinkel and Rao [8], Inspired by a 
work of Edmonds and Fulkerson [51 on general bottleneck problems, they called the procedure 
the Threshold Method. The algorithm can be described in the following way: First determine a 
"good" lower bound z for the optimal objective value z. Then define a network rfiz) in the fol¬ 
lowing way: Let V - (s,/}0 WOJV with M — [1,2. m } and N » {m + 1. m + n}. 

be the nodeset of X(z). Now every node / € M is connected with s by an arc, respectively, 
each m + j € N with t. A pair (i,m + j) is connected by an arc only if t :j < z holds. 

The arc-capacities d k i are defined by 

dyi ai U — 1 . m) 

d ,.m+j °° 0 - 1. m\j - 1. n) 

d m +j.i '■= bj 0-1. n). 

Now a maximum flow / - (f 0 ) from source s to sink t is determined using the labeling method 
of Ford and Fulkerson [6], If for the maximal flow value v — £ a, holds, the optimal solution 

for the BTP is obtained defining 

x ij /i.m+y O' - 1. m\j - 1, .... n) 

If v < I a, the labeling procedure yields simultaneously a minima! cut (X,X) in S(z) 
with s € X and / € X. This cut has the property 

t 0 < r ( t,m + j) V ( X,X) 

To obtain a solution to (1.1M1.3) it is therefore necessary to use an arc ( i,m + j) with t tJ > z. 
Hence, we determine 

z* :— min {r /y |/ € A', m + j €X) > z- 
Define z :- z* and repeat the process. 

After (m • n) iterations, at most, an optimal solution is obtained. The algorithm is sum¬ 
marized in the following flow-chart: 
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Figure 1. Flow-chart of Ihe Hungarian method 


During the algorithm the network X(z) need not be constructed explicitly. All operations of 
the labeling procedure can be performed directly on the "time" — matrix T - (t 0 ). 

3. A PRIMAL METHOD 

The first method for solving BTP proposed by Barsov [l] is a primal algorithm which 
proceeds "dual" to the Hungarian method in some sense. In the course of the Hungarian 
method the "threshold" z is successively increased until a feasible solution can be found. In this 
method a feasible solution is always at hand and the threshold-value is successively decreased 
until no further solution can be found. Starting from a feasible solution x - (x 0 ) with objec¬ 
tive value z a cost-matrix D » (dy) is defined by 

| 0 if l u < * 

d'J " { 1 else 

Now a TP with cost-matrix D is solved. If the optimal solution y *» Oty) has an optimal 
value z(y) > 0 the solution x — (x y ) is optimal for BTP, Otherwise 

z* max { t,j ly^ > 0} < z holds. 


Define z z* and x y and repeat the process. 

Computational experience shows (8) that this method is inferior to the Hungarian 
method. Recently, Finke and Smith 17] developed an improved primal method. 
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Starting from a lower bound r and a "good" start solution x - (x„), a cost-matrix 
D - (4,) is defined by 


2 + 1 


with lx] greatest integer < x. 


Now the associated TP is solved using the well known MODI-method and x — (x</) as 
start solution. If the optimal solution y — 0^) has objective value z{y) — 0, then y is optimal 

for BTP. Otherwise, we consider the dual variables i — 1, .... m and v y , j — 1. n 

associated with the optimal solution y. Then 

u, + v , ^ d,j i = 1, ... , m \ j — 1. n 

y.j > 0 -*• U, + Vj - d,j 
and 

z 'min lt 0 \u, + \j > 0} > z holds. 


It can be shown that z’ is a lower bound for the optimal value for the BTP. Define z :« z' 
and repeat the process. In an implementation the actual basic solution of the TP is stored as a 
tree using the list structures proposed by Glover and Klingman [91, Srinivasan and Thompson 
U2). Due to this technique primal methods are superior to other methods in the case of TP. 

The following flow-chart summarizes the algorithm of Finke and Smith. 



Figure 2. Flow-chart or the primal method 
of Finke and Smith 
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4. THE SHORTEST AUGMENTING PATH METHOD 

This method was proposed by Tomizawa 114} for solving TP and applied to BTP by Derigs 
and Zimmermann [4], A theoretical foundation of this method can be found in [2J. During 
the procedure so-called /c-subproblems (P k ) have to be solved for 1 < k < m. 

(P k ) min max (///lx*/ > 0} 

„ zz, 1=1,2. k 

£ X 'J “ 0 else 
y-i 

£ Xu < bj j - 1 . n 

I-1 

Xy > 0 i — 1. m \ j — 1. n. 

The solution for (P x ) is obvious. Starting from an optimal solution x - (x„) for ( P k ) 
with k — 1, , m - 1 an optimal solution for (P* +1 ) is obtained by means of augmenting 

paths. For this purpose we partition (1,2. n] — A/, 0 jV 2 with 

A, — |y | x,j - byj "saturated columns." 

Now a sequence Ax of mutually distinct matrix-entries 

(* + W. UlJ\), (/iJj), U 2 .j 2 ), .... U,J r ), (Wr+I> 
is called an augmenting path with respect to x - (x y ) if 

x. , > 0 for q - 1.2. r and 

y, € N, for q - 1.2, ... ,r and j r +\ € N 2 . 

The capacity of Ax is defined by 

k 

cap (Ax) min {a* + ,, b Jr+{ - £ x,. >f+| , «} with 
« min {x, . 14 - 1,2 .r}. 

9 ,J 9 

Let us define x x © Ax by 

Xjj - cap (Ax) for (/.y) - (/|J|). U,,j r ) 

x,j Xy + cap (Ax) for (ij) € Ax\{(i,,y,)l 9 - 1, .... r\ 

Xy else 

and 

w (Ax) max [t,j\x tJ > 0). 

Let D k be the set of all augmenting paths with respect to x - (x u ) and Ax € D k with 

(4.1) w(Ax) < w(Ax) for Ax € D k 

(4.2) cap (Ax) - a* + j 

then x x © Ax is an optimal solution for (P k . M ). 
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An augmenting path Ax € D k with property (4.1) can be determined using a modified 
Dijkstra method. If Ax does not fulfill (4.2) we split (P k +\) into two subproblems (P k +t ) and 
(P k + 1 ) with <j* + i cap (Ax) and a k+l a k+) - a k+l . Thus, (P k +i) can be solved in a finite 
number of steps. 


Computational experience shows that the method can be improved with the aid of the fol¬ 
lowing starting procedure. First we determine a "good" lower bound z for the optimal value z 
and a "good" partial solution x - (Xy), i.e. 


(4.3) £•*//< a, , ' - 1. m 

(4.4) Jjxy < bj ,7-1 . n 

(4.5) Xy > 0 , i - 1,-w, 7 - 1 . * 

with the property 


(4.6) Xy > 0 -► fy < * 

such that 

£ £ Xy is large. 

/-i 7-1 

A partial solution x with (4.3) • (4.6) is called ^-feasible. Computational tests have shown 
that it is not useful here to determine a z-feasible solution with maximal possible sum. Some 
simple heuristics will give solutions which are within 5% from optimality much faster and, with 
respect to the running time of the complete procedure, this is of advantage. The complete pro¬ 
cedure is summarized in the following flow chart. 



Fioline 3. Flow-chart of the shortest 
augmenting path method. 
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5. LOWER BOUNDS AND STARTING PROCEDURES 

All three methods presented above start with a lower bound z for the optimal value z. 
Computational experience shows that the bound proposed by Garfinke) and Rao is favorable. 
Row-thesholds z, and column-thresholds z j are computed with 

z, min max {|jc 0 > 0} subject to 

I *u > a , 

j -1 

0 ^ X/j < bj for / = 1. m 

and z J defined analogously. 

Then z max {zi ,z 2 . ••• . z m , z',z 2 , .... z") is a lower bound. This bound is easy to 
calculate and yields a good estimation for the optimal value in most cases. 


The primal method of Finke and Smith heeds a good starting solution x = (xy). They 
propose the following method. First a z-feasible partial solution x — (xy) is constructed. Such 
a partial solution is used for the shortest augmenting path method, too. A z-feasible solution 
can be obtained in the following way. For every row / (column j) the number 7] ( P) of feasi¬ 
ble matrix entries (i,j) subject to t tJ < z is computed. Now the row i 0 (column j Q ) with 
minimal 7) # > 0 (7 J ° > 0) and positive a /g (by g ) is determined. If no such row (column) exists 
the partial solution is completed. Otherwise, determine in row /„ (column j 0 ) the feasible entry 
('Wo) with minimal number 7 J ° (7^) and positive b Jg (a ig ). Then define 


min \a,. bj 



a '0 

% ~ x 'o-jo‘ 


Update T : and P and repeat the process. 


This way a good z-feasible partial solution x - (x tJ ) is determined. If x — (xy) fulfills 

0.1) and (1.2) then it is an optimal solution for BTP. If the partial solution is not optimal in 
the second step of the starting procedure proposed by Finke and Smith, the remaining supply is 
distributed using a north-west-corner rule. Finke and Smith call their procedure threshold- 
totals-method. 


The following flow-chart summarizes the starting procedures for the three methods. 
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Figure‘ 4. Flow-chart of the starting procedures. 


6. COMPUTATIONAL RESULTS 

FORTRAN IV implementations of all three methods were tested on a CDC CYBER 76 of 
the Computer Center of the University of Cologne. We used the following codes: 

HUNGAR - a version of the Hungarian method which Professor Garfinkel kindly made 
available to us 

PRIMAL - a version of the primal method which Professor Finke kindly made available 
to us 

SAP - an improved version of the BTP-code which is listed in [4], The improve¬ 
ment concerns mainly the starting procedure. 

The HUNGAR-code uses the maximal cardinality NL of expected postive x u 's in any 
column as input data. The running time of the program is highly dependent on the choice of 
NL since the value and position of the x,/s are stored in NL x N matrices. Thus, the choice of 
smaller "NL" decreases the running time. But if the actual number exceeds NL, termination 
occurs. Choosing NL " M/2 all problems were solved. 

We considered rectangular as well as quadratic examples. The integer coefficients of the 
time matrix 7", the supply vector a and the demand vector b were generated by a machine 
independent uniformly distributed pseudo random number generator in the interval [1,2 31 -1]. 
Afterwards these numbers were transformed into the intervals [1,6] with 
b - 10, 100, 1000, 10000, 2 31 -1. Twenty-five examples were generated for each combina¬ 
tion to calculate the mean running time. Numerical experience showed that the running time 
of HUNGAR and SAP for solving rectangular problems with m > n is less than for the 
equivalent transposed problem with m < n. Therefore, we recommend that rectangular prob¬ 
lems always be solved in the form with m > n. 
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The following table shows the mean running time for quadratic problems with 
m — n — 100. 


TABLE l — Mean Running Time in CPU-Seconds for (100 x 100) BTP 


a„bj 


1-10 

1-100 

1-1000 

1-10000 

l-<2 31 -1) 


SAP 

.112 

.177 

.177 

.192 

.181 

2 

PRIMAL 

.106 

.124 

.124 

.133 

.126 


HUNGAR 

.226 

.400 

.383 

.478 

.504 

8 

SAP 

.107 

.120 

.144 

.171 

.150 


PRIMAL 

.094 

.106 

.105 

.106 

111 


HUNGAR 

.332 

.582 

.611 

.588 

.672 

© 

SAP 

.089 

.140 

.188 

.121 

.223 

2 

PRIMAL 

.092 

.098 

.119 

.106 

.114 

- 

HUNGAR 

.333 

.439 

.548 

.603 

.713 


The mean running time of PRIMAL is significantly better than the other's. For SAP this 
is simply caused by a single "ill conditioned" problem for which the running time is more than 
tenfold the running time of the other 24 examples of the group whose running time is compar¬ 
able to those of PRIMAL. 

SAP and PRIMAL are shown to be relatively insensitive to the range of the parameters 
Oj.bj and t 0 while the running time of HUNGAR is more than doubled when the range for r„ is 
increased from b — 10 to * - 2 31 -1. 

Then we modified the problems subject to 

/1 j — /,) — 0 for / “ 1. m and j — 1. n and 

a\ - f>t - max {max{a,}, max (&,}}. 

For these perturbated problems the starting procedure has no effect and thus the computational 
behavior of the "pure" algorithms is shown. 

Table 2 shows that PRIMAL and HUNGAR are more dependent on the quality of the 
starting procedure and the range of the time values t,j than SAP. 


TABLE 2 — Mean Running Time in CPU-Seconds for the Perturbated 
(40 x 40) Problems 


Oi.bj 


1-10 

1-100 

1-1000 

1-10000 

l-(2 31 -1) 

© 

SAP 

.014 

.105 

.122 

.117 

.121 


PRIMAL 

.015 

.118 

.369 

.446 

.448 


HUNGAR 

.016 

.269 

.608 

.720 

.686 


SAP 

.014 

.153 

.196 

.202 

.169 

$ 

PRIMAL 

.014 

.133 

.389 

.483 

.504 

4 

HUNGAR 

.026 

.436 

.790 

.922 

.988 

g 

SAP 

.014 

.161 

.179 

.217 

.164 

© 

PRIMAL 

.015 

.132 

.398 

.483 

.506 

- 

HUNGAR 

.025 

.422 

.834 

.937 

.983 
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Then we applied the algorithms to bottleneck assignment problems i.e., problems with 
a, — b, = 1. 


TABLE 3 — Mean Running Time in CPU-Secondsfor (100 x 100) 
Bottleneck Assignment Problems 



'U 

1-10 

1-100 

1-1000 

1-10000 

1 - (2 31 -1) 

| E 

SAP 

.074 

.074 

.074 

.073 

.071 

E u 

c X) 

00 o 

PRIMAL 

.083 

.092 

.100 

.096 

.105 

tfl >* 

CO CL 
< 

HUNGAR 

.055 

.103 

.143 

.140 

.155 


For assignment problems SAP needs at most n/2 iterations while PRIMAL will perform 
a number of degenerate pivot operations. Again PRIMAL and SAP show a more robust nature 
with respect to the range of the time parameters t,j than HUNGAR. 

More computational results can be found in [3] and [4]. 

7. CONCLUSIVE REMARKS 

In this paper we presented three basic methods for solving the bottleneck transportation 
problem. It is beyond the scope of this study to give an entire review of all the different 
methods which are available for solving BTP. Nevertheless, we think that the presented pro¬ 
cedures are representative in the sense that the main approaches for tackling combinatorial 
optimization problems are specialized to the bottleneck transportation problem. 

The computational tests indicate that PRIMAL and SAP outperform HUNGAR in gen¬ 
eral. The optimal choice between PRIMAL and SAP is data dependent and should be specially 
made for every problem to be solved. If it can be expected that the start-heuristic is perform¬ 
ing well, PRIMAL is recommended. Otherwise, SAP is of advantage. 

We want to close our discussion by referring to another computational study. Werner 
(15] compared versions of SAP and PRIMAL on special problems the data of which was 
defined by some structured transportation networks incorporating time tables for transportation 
vehicles, waiting-times and different regions with different magnitudes of demand and supply. 
Werner reports that the running time for PRIMAL was highly dependent on the distribution of 
data and the more the data was not uniform the more SAP became superior. 
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SEQUENCING n JOBS ON TWO MACHINES WITH SETUP, 
PROCESSING AND REMOVAL TIMES SEPARATED 
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ABSTRACT 

The paper extends the machine flow-shop scheduling problem by separating 
processing time into setup, processing and removal times. 


In a recent paper Yoshida and Hitomi [2] extended the classic two machine scheduling 
problem first introduced by Johnson (1 ]. They allowed setup times to be independent of the 
processing times. This paper is the further extension of the Yoshida and Hitomi model; it 
allows for separation of processing time into setup time, processing time and removal time for 
each job on each machine. 

For example, consider a machine shop environment. Operations associated with each job 
on each machine when the machine/operator is available could be summarized as follows: 

1. Setup time independent of the unit to be processed. This operation consists of activi¬ 
ties such as obtaining the blueprints, procuring the necessary tools, fetching the 
required jigs and fixtures and setting them on the machine. 

2. Setup time that is a unit dependent. This operation includes the time required to set 
the unit in the jigs and fixtures and to adjust the tools as required. 

3. Processing time. 

4. Removal time dependent on the unit. This operation includes the times for activities 
such as disengaging the tools from the unit, and releasing the unit from the jigs and 
fixtures. 

5. Removal time independent of the unit. This operation includes activities such as dis¬ 
mantling the jigs, the fixtures and/or tools, inspecting/sharpening of the tools, 
returning them to the central depository, cleaning the machine and the adjacent area. 

Since activities 2, 3 and 4 are unit dependent, their times could be combined and desig¬ 
nated as the processing time. However, times for the activities I and 5 are independent of the 
unit to be processed. 
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Formulation of the Problem 


Let, 


P" - processing time of job ion machine m. This includes times for the activities 2, 3 and 
4 described before, i - 1,2, ... n, m - 1,2. 

Sr - setup time independent of the unit, i.e., activity 1, for job / on machine m. 

R” —removal time, i.e., activity 5, for job / on machine m. 

C" —completion time of job / on machine m. 

In addition, in the mathematical development, it is assumed that the jobs are designated 
so that the job in the Ah position of the processing sequence is job /. Since the setup is 
independent of the unit, if there exists an idle time on machine II, setup on machine II can be 
done before the unit is available from machine I. The unit is available from machine I when 
processing on the unit is completed; however, the machine is not available for the next job 
until the removal operation on the machine is finished. Figure 1 shows these activities graphi¬ 
cally. 




— a VMM 


Figure l. Graphical illustration of two machine problem 


Completion time for job i on machine I is given by, 

0) c,' - s' + p j + 

Completion time for job »' on machine II is 

(2) C, 2 - (C,L, + 5/ + P,') + P, 2 + Rj 2 

if C 2 - 1 + S 2 < C ( L| + S' + Pf (see Figure 1) 

or 

C, 2 - C ( i, + S, 2 + P} + R 2 
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C,l, + 5, 2 > C,L, + 5,' + />'. 


The combination of (2) and (3) gives 
(4) C, 2 = Max (C, 1 - R,' - S, 2 . C,l, ] + S 2 + P 2 + R, 2 . 

By successive application of (4) using (1), the total elapsed time is given by 


(5) 


T = C 2 = Max 


Max 

O^u^n 


£ (5/ - S 2 + /»/) - £ (/>, 2 + 2?, 2 - /?,') 


, 0 


+ T (S 2 + P 2 + R 2 ). 

The optimal schedule can be found by using Johnson’s method to solve a two machine flow- 
shop problem where the processing times of job ion machines I and II are S, 1 — 5, 2 + P and 
P 2 + R, 2 — R t \ respectively. 


I 
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A NOTE ON THE MAXIMIN VALUE OF TWO-PERSON 
ZERO-SUM GAMES* 


Alvin E. Roth 


Mellon Professor of Economics 
University of Pittsburgh 
Pittsburgh, Pennsylvania 


1. INTRODUCTION 

The theory of two-person, zero-sum games has occupied a special place in the general 
theory of games ever since its introduction by von Neumann and Morgenstern [17]. They 
presented the theory of two-person, zero-sum games as a natural extension of the theory of 
rational choice under uncertainty, as modeled by the maximization of an expected utility func¬ 
tion. Their conclusion was that rational players should always play their maximin strategies in 
such games, and should regard the maximin payoff as the "value" ol the garne t 

Subsequent authors have questioned both the validity and generality of these conclusions. 
Ellsberg [4], for instance, points to gaps in the arguments which von Neumann and Morgen¬ 
stern present in support of their conclusions, and raises questions related to the interpretation 
of the naure of the outcomes of the game, while Aumann and Maschler [1] discuss issues 
which arise in passing from the extensive to the normal form of a game. This latter discussion 
has generated a lively controversy (cf. Aumann and Maschler [21, Davis [31, Owen [71, Taylor 
[131). McClennen [6] also examines gaps in the arguments presented by von Neumann and 
Morgenstern, and, like Aumann and Maschler, concludes that the prescription that rational 
players should choose maximin strategies cannot be derived directly from the principle that 
rational individuals are utility maximizers. 

There is thus a gap between the rationality assumptions which insure that an individual 
evaluates single person decision problems as a utility-maximizer, and the assumptions which 
insure that he evaluates two-person, zero-sum games in terms of the maximin value. Two dis¬ 
tinct issues are involved here, since the maximin value is intended to be used both to evaluate 
alternative games and to identify maximin strategies as "rational" choices. This paper addresses 
the first of these issues. 

One approach to studying the maximin value is the axiomatic approach of Vilkas [151 and 
Tijs (141, who consider arbitrary functions defined on games and present axioms which are 


•This work was supported by National Science Foundation Grant SOC75-21820 to the Institute for Mathematical Stu¬ 
dies in the Social Sciences, Stanford University, and by National Science Foundation Grant SOC78-09928 to the 
University of Illinois. 

tStrictly speaking, this conclusion applies only to games having saddlepoints in the set of admissable strategies. Such 
games are sometimes referred to as determined games. 
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uniquely satisfied by the maximin value. This paper takes a different, but related, approach, 
and shows that the maximin value of a two-person, zero-sum game can be interpreted as a 
rational individual’s von Neumann-Morgenstern utility for playing the game, and demonstrates 
necessary and sufficient conditions on the individual’s preferences over games for this to be the 
case. That is, we will consider conditions which an individuals's preferences must obey in addi¬ 
tion to those which insure that he is a utility maximizer, in order for his utility function to coin¬ 
cide with the maximin value. This approach has been used to study other classes of games: 
e g., cooperative games with side payments (Roth (81, (9]), simple games (Roth [10]), and bar¬ 
gaining games without side payments (Roth [11], [12]). The concluding section of this paper 
will discuss briefly the relationship between the results obtained here and these other results. 

It should perhaps be emphasized at this point that it is not the purpose of this paper to 
further explore the mathematical properties of the maximin value, which are well understood. 
The questions which have been raised in some of the papers cited above concern, rather, the 
interpretation which can be placed on the maximin value as part of a theory of rational (i.e., 
utility-maximizing) behavior. The purpose of this paper is to further explore this latter issue. 

2. THE MAXIMIN VALUE AS A VON NEUMANN-MORGENSTERN UTILITY 

A two-person, zero-sum game g will be denoted by a triple g = (S u S 2 ,u), where S| and 
S 2 are arbitrary sets, and u is a function such that u:S | x — R 2 , and such that, for any 
(s,f) € S, x S 2 , u,(s,r) = -u 2 (s,t). The interpretation is that in the game g , the player in 
position / (/ =1,2) will have available a set of strategies* S„ and will compete for prizes over 
which his preferences are represented by the utility function u,. His opponent’s preferences 
will be precisely opposite his own. (Note that in general it is not necessary that the players 
compete for a single set of physical prizes.t It is sufficient that the results of the game be such 
that the two players are never in agreement over which of two outcomes is preferable.) 

This formulation defines the game independently of the individuals who will play it. Con¬ 
sequently, it is meaningful to consider the preferences of some individual over the positions in 
a game, or in different games. Letting G be some class of two-person, zero-sum games, and 
letting N = {1,2) be the set of positions, we will be considering an individual’s binary prefer¬ 
ence relation P defined on IV x C, the set of positions in a game. For i, j € N and g , g' € G , 
the statement ( i.g)P(j,g') should be read "it is preferable to play position / in game g than to 
play position j in game g'." We will interpret P as a strict preference relation, and define an 
indifference relation I and a weak preference relation W in the usual way.* We will assume that 
the preference relation is defined over the mixture set of lotteries generated by N x G, and 
that it satisfies the rationality conditions necessary for the existence of an expected utility func¬ 
tion v. 


’Typically the strategy set S, could be the set of all probability mixtures of some finite set of pure strategies 
+Two person, zero-sum games are sometimes interpreted as games in which, for every pair of strategy choices, one 
player’s gains are precisely equal to the other's losses in terms of some physical commodity, such as money Under 
this interpretation, the game would only be zero-sum in the rare case that the two players have precisely opposed 
preferences for lotteries involving money In our formulation, the game specifies the utility payoffs to be faced by the 
players For a given pair of individuals to play a specified game, the monetary rewards to each would have to be ad¬ 
justed in terms of each individual's utility function for money 

*So, for a.b € N x G. alb if and only if neither aPb nor bPa. and if and only if either aPb or alB. 
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Since we are interested in preferences over a possibly infinite set of alternatives (e.g., let 
G be the set of all matrix games) it will be simplest to refer to the axiomatization of utility due 
to Herstein and Milnor IS], rather than to von Neumann and Morgenstern’s original axiomati¬ 
zation. Denote by [p(i,g)\ (1 — p)(J,h)] the lottery which with probability p has an individual 
play position i in game g, and with probability 1 - p play position j in game h. Then we assume 
that the preference relation is defined on a mixture set containing all such lotteries, and that 
this preference relation has the properties of transitivity, continuity, and substitutability 
sufficient for the existence of an expected utility function (Herstein and Milnor [5]). That is, 
there exists a function v (unique up to origin and scale) which preserves preference (v(/,g) > 
vO',/i) iff ( i,g)P(j-h )) and evaluates the utility of a lottery by its expected utility (v[p(/,g); 
(1 - p)(J.h)] = pv((/,g)) + (1 - p)\((j,h)). 

The fact that we are considering preferences defined over games, which are themselves 
defined in terms of an expected utility function, puts some additional restrictions on the allow¬ 
able preferences. In particular, for any real number k , let g k be the (degenerate) game defined 
by g k -* (Si,S 2 ,u) such that |S)| - |S 2 I “ 1. and «|(S],S 2 ) - k. Thus g k is the game* in which 
each player has only one strategy to choose, which results in the player in position 1 receiving a 
utility of k , and the player in position 2 receiving a utility of — k. Since k is an expected utility, 
we know precisely how an individual’s preferences should behave over alternative games g k . 
r Atmally, we assume that the set G contains the set of degenerate games, and impose the fol- 
. /ing conditions on preferences over degenerate games, which have the effect of embedding 
the space of utilities into the space of games. 

(1) (\. gl )P(l.go). 

This simply states that it is preferable to play position I in the degenerate game g| which 
awards the player in position 1 a utility of 1 than to play position 1 in the game g 0 and get a 
utility of 0. 

(2) For any real numbers cand A: such that c ^ 1, 



The first part of this condition states that a player is indifferent between receiving a utility of k 
for certain in the degenerate game g k , or participating in a lottery which gives him a utility of ck 
with probability 1/c and a utility of 0 otherwise. The second part of the condition states that a 
player is indifferent between receiving a utility of k or —fc, each with probability 1/2, or receiv¬ 
ing a utility of 0 for certain. This is precisely what is meant by the statement that the function 
u is an expected utility function. 

Condition (1) permits us to normalizet the utility function v in the natural way, and we 
henceforth take v(l,g 0 ) - 0, and v(l,g|) - 1. Condition (2) then completes the embedding of 
utilities into the space of games, giving us the following result. 

PROPOSITION 1: For any real number fc, v(l,g*) — k. 


•Technically speaking, our definition allows there to be more than one game g k , since we could allow different one- 
element strategy sets. This difference is inessential, and can be ignored. 

(Recall that v is determined only up to choice of origin and scale. 
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PROOF: Since v is a utility function, the first part of (2) implies that, for c ^ 1, 
- vll/c(l,g r *); (1 - (l/c))(l,g 0 )]. But v is an expected utility function, so the utility 
of a lottery is its expected utility; i e., v{l/c(l,&*); (1 - (l/c))(l,£o)l “ d/c)v(l,g ck ) + 
(1 - (l/c))v(l,£ 0 ) ” (l/e)v(l,g rt ). Thus, for any k , v(l,g rt ) - cv(l,g t ) for c > 1. For 
0 < c < 1, let d - 1/c, and k' - ck. Then vO.g*-) « dv(\,g k -), so vU,&*) - cv(l.g*) for 
any k and any c > 0. In particular, for any k > 0, vO.g*) — fcv(l,£i) — k. But the second 
part of (2) implies that v(l,&) - -v(l,£_*) for any k , and so v(l,£*) - k for all k , as was to 
be proved. 

So far we have imposed conditions on preferences over degenerate games which reflect 
the fact that they are associated with the underlying expected utility function, but we have yet 
to impose any conditions which reflect the zero-sum nature of all games in G. The following 
conditions will suffice for our purposes. 

(3) For all g € G, (l,g 0 ) 

This condition states that a player is indifferent between getting a utility of zero for certain in 
the degenerate game # 0 or participating in a lottery which gives him an equal probability of 
playing either position in any two-person, zero-sum game g. The motivation for this require¬ 
ment is the idea that if a given zero-sum game g yields an advantage to a player in one of the 
positions, then it must yield a corresponding disadvantage to the player in the other position. 
An immediate consequence of condition (3) is the following: 

PROPOSITION 2: For all g € C, v(U) - -v(2 ,g). 

PROOF. (3) implies that v(U 0 ) - vll/2(l,*); l/2(2,g)] - (l/2)v(l,s) + (l/2)v(2,g). 
But v(l.s 0 ) - 0, so v(l.g) - -v(2.#). 

To state the next condition, consider a game g - ( S\,S 2 ,u ) and a given position i € N. 
for every strategy s € S,, define k(s) - kig.i.s ) - inf u,(s,r), where j ^ /'.* Then in addition 

to the conditions already imposed on the preferences, we require the following. 

(4) For any g € G, i € N, and s € S„ ( i,g ) W(i,g k{s) ). 

If the game g were such that S, contained only a single strategy, then condition (4) would sim¬ 
ply be a version of the "sure thing” principle, which states that any prospect is at least as desir¬ 
able as the worst outcome which can result from that prospect. As it stands, the condition also 
reflects that individual i is free to choose any strategy s in S,. Condition (4) implies the follow¬ 
ing proposition. 

PROPOSITION 3: For any g € G, 
v(l.g) ^ sup inf u,(s,(). 

3(S, i€S 2 



•So kif?. l.s) - inf M.(j.r), and k(g, 2,s) - inf 
t(S 2 r(S j 


y (1,*); y a*) 
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That is, the utility of playing a game is at least as large as the "security level" which the game 
provides. No further assumptions are needed to permit us to reach the following conclusions. 

PROPOSITION 4: For any g € G, 

inf sup «|(s,r) > v(l,g) > sup inf uAs.t). 
i(S 2 s«s, jeS, its 2 

PROOF: The right hand inequality is simply the first part of Proposition 3, while the left 
hand inequality follows from Proposition 2 and the second part of Proposition 3, together with 
the fact that u,(s,f) - -u 2 (s,t). 

A game g is said to have a saddlepoint if there exist strategies s € S, and t € S 2 which 
achieve the inf sup and sup inf of Proposition 4, and make them equal. For such games. Pro¬ 
position 4 has the following immediate consequence. 

PROPOSITION 5: If g € G has a saddlepoint, then 

v(l,g) - max min u.(s.r) — min max ui(s,/). 

s€S, tiS 2 its j jES, 

The famous Minimax Theorem of von Neumann [16] establishes that if g is the mixed exten¬ 
sion of a finite two-person, zero-sum matrix game, then g has a saddlepoint and Proposition 5 
applies. 

DISCUSSION 

The previous section establishes conditions on an individual's preferences over games 
which yield the maximin value as his utility for a given game. Insofar as these conditions are 
consistent with our other notions about zero-sum games, this shows that the maximin value is 
consistent, both formally and substantively,* with the notion of rationality embodied in utility 
maximization. The "machinery" necessary to establish this result shoiuld not, however, be per¬ 
mitted to obscure the essential simplicity of the argument. 

The two substantive conditions on preferences are (3) and (4). Condition (4) serves to 
establish a lower bound on the utility of playing the game from either position. Condition (3), 
by fixing the utility for one position in the game with respect to the utility for the other posi¬ 
tion, combines with (4) to produce an upper bound on the utility. For a game with a 
saddlepoint these bounds coincide, and so the utility is determined. 

Since the maximin value is often interpreted as being an unnecessarily conservative 
assessment of a game (cf. Ellsberg’s [4] "reluctant duelist"), it is interesting to note that condi¬ 
tions (3) and (4) actually prevent the utility function from reflecting an excessively conserva¬ 
tive attitude. In particular, if g is a game without a saddlepoint, then the utility of playing g in at 
least one of the positions must be strictly greater than the maximin value for that position. 
This is because Proposition 2 requires that v(\,g) - —v(2 ,g) y and, in a game without a 
saddlepoint, the maximin values for each position are not the negatives of one another. 


‘Any numerical index Tor certain events can of course be extended to a function which preserves expected value over 
lotteries. The purpose of this paper is to investigate what substantive assumptions about preferences need to be made 
in order to interpret this extension as an expected utility function. 
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Note that if the class G does include games which have no saddiepoint, then the condi¬ 
tions imposed in the previous section do not completely determine the utility function v. To 
do so, it would be necessary to precisely specify a player's attitude towards the additional uncer¬ 
tainty present in games with no saddiepoint. The same phenomena occurs in studying coopera¬ 
tive games, in which it is necessary to specify a player's attitude towards strategic risk (cf. Roth 
191). 


As a Anal observation, note that, although the maximin value is linear in degenerate 
games, it is not linear in arbitrary games. Let g and h be finite matrix games of the same 
dimensions, let p be a probability (i.e., p € (0,1)), and let g' — pg + (1 — p)h (where we 
employ the usual conventions of matrix arithmetic). Then each element of the matrix g' is the 
expected value of the corresponding element of the random matrix defined by the lottery 
[pg\(\ - p)h]. We might cot\jecture that a player would be indifferent between playing a given 
position in the game g' or taking the same position in the lottery. But for preferences obeying 
the conditions considered in this paper, this is not the case. That is, v[p(r,g); (1 - p)(i,/i)l 
equals pv(i,g) + (1 — p)v(i,h) and is in general not equal to v(i,g'), for i — 1,2.* In this 
respect, the preferences considered here over zero-sum games resemble the preferences over 
bargaining games considered in Roth [111, [121 more than the preferences over games with 
sidepayments considered in Roth [8], (9). 
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ABSTRACT 1 

Certain types of communication nodes can be viewed as multichannel ! 

queueing systems with two types of arrival streams Data arrivals are character¬ 
ized by high arrival and service rates and have the ability to queue if all service 
channels are busy Voice arrivals have small arrival and service rates and do 

not have the ability to wait when the channels are full. Computational pro- 1 

cedures are presented for obtaining the invariant probabilities associated with 
the queueing model 


1. INTRODUCTION 

The steady growth of the data processing industry and the telephone network over the last 
Tour decades has introduced computer-communications networks as efficient transportation 
vehicles for the remote sharing of information data bases. Recent Defense Communications 
Agency studies have shown the desirability of a network which integrates voice, interactive, and 
bulk data for the 1980's (Rosner [10], Schmitz, Saxton, Huang and White [II]). Several addi¬ 
tional studies relating to information processing growth in the next few decades portend new 
data/voice services with substantially increased data flows (Rich and Schwartz [8], Rosner [9]). 
There are several different time division integration methodologies proposed for combining 
voice and data demands over the available channel bandwidth (Coviello and Vena [2], Fischer 
and Harris [4], Frank and Gitman [5]). A competitive allocation scheme allows the data and 
voice calls to compete for time slots using a first-come, first-served scheme. 

Simulation models have been developed to study competitive allocation schemes; how¬ 
ever, because of their expense and the fact that rare events are of interest, simulation studies 
are not necessarily appropriate. Analytical models have also been developed to describe com¬ 
petitive allocation (Bhat and Fischer [1], Fischer [3], Fischer and Harris [4]). The analytical 
models developed to date all share the common drawback of not being tractable for problems of 
reasonable size. The purpose of this paper is to present an analytical model of the dynamic 
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movement of voice and data calls within an integrating device and give a tractable solution pro¬ 
cedure for the model. The allocation scheme is that used by Bhat and Fischer [1], The 
methodology for the solution scheme combines procedures developed by Neuts [6,71 and 
Wong, Griffin, and Disney (121. 

2. THE MODEL 

There are two classes of calls that arrive to the service center: data calls and voice calls. 
The arrival process of data calls is assumed to be Poisson with rate X) and the arrival process of 
the voice calls is assumed Poisson with rate X 2 ; the two arrival processes being independent. 
Although each data call actually consists of a random number of discrete bits, the length of 
each data call is assumed to be adequately approximated by a continuous exponential random 
variable with mean l//xi and the length of the voice call is assumed exponential with mean 
1 /fi 2 - The characteristics of data calls are frequent arrivals and short length whereas the voice 
calls have infrequent arrivals with a long service time. A buffer is available that can store 
incoming data calls when all channels are busy. The size of the buffer is such that no effective 
limit is placed on the queue size. Voice calls, however, may not form a queue so that any 
incoming voice call is lost whenever all channels are busy. The service center has c channels 
using a FIFO discipline. 

The data/voice integrating system is modeled as a Markov process with a two dimensional 
state space given by 

E — {(n.m): m = 0,1.c and n = 0,1, ... 

The vector p (ordered lexicographically) denotes the steady state probabilities where p(n.m) is 
the steady state probability that there are n data calls in the system and m voice calls in the sys¬ 
tem. For ease of notation the vector p is partitioned as 

(1) P - (Po. Pi, ■ ■ • ) 

where the mth component of p„ is p(n,m). The vector p is the solution to the equation 

pQ - o 

(2) pl-1 

where 1 and 0 are vectors of all ones and zeros, respectively, and Q is the transition rate matrix 
for the Markov process representing the data/voice system. The matrix Q is infinite and is 
given by 


A 0 A 
M | A, A 


o 


(3) 


Q - 


M c _ 


i 


o 


A c -1 A 

Me— i Aq A 

M c A c A 
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where the submatrices are defined, for i.j.n — 0,1. c, by 

lx, if 1-7. 

^ A(/J) — | q otherwise; 

nfi | if i - j and j < c — n, 

(5) - (c — j)fi ] if i — j and ./ > c — n, 

0 otherwise. 

k 2 if i — j - 1 and i < c — » — 1, 

Hi if i - j + 1. 

( 6 ) if/ , 7 

0 otherwise; 
where a„(i) is such that the row sum is zero. 

3. SOLUTION PROCEDURE 

The general structure of the matrix Q is identical to the type investigated by Neuts [6] 
and similar to the finite system investigated by Wong, Griffin and Disney [12]. The solution 
method as given by Neuts [6] is a two step process. First, it is necessary to find the matrix R 
such that 

(7) R 2 M c + RA c + A - 0. 

The matrix R gives the relationship 

(8) p c+ * - p f _,«* +l for X = 0,1. 

The second step is to let p = (p 0 .P f -i) and then solve for p u'i^g equations 

(9) p T - 0 
and 

p 1 + p f _, R (/ - R)-' I - 1 


M c -\ “F R M c 


It will turn out that for this communication problem, the matrix R of equation (7) is easy to 
obtain whereas jk from equation (9) gives computational difficulties which will be overcome by 
using a technique utilized by Wong, Griffin and Disney [12]. 
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The matrix R is lower triangular since the matrices M c , A c , and A are lower triangular. 
Therefore, the diagonal elements of R 2 are simply the square of the diagonal elements of R 
Thus, the following equation results when considering the diagonal elements in (7) for 
/ - 0.c 

(10) r (c — i)fi, — /•„[*] + ifi 2 + (c — /)/i|l + A| = 0. 

The solution of (10) yields 

r cl — A |/(cm 2 + A|) 
and for / = 0.c — 1 

r H — {A | + /'m 2 + (c — /)mi — ((A i + /'M 2 + (c — »')mi) 2 

(11) -4A,(c - /')mi1 I/2 )/|2(c - »)mi). 

The oflf-diagnonal elements are not quite as straight forward but by considering the ij element 
with i > j, equation (7) yields 

<c - ./>Mi L r, k r kl + ]£ r lk A t (k,j) - 0 

A-y *-0 

which in turn yields 

/-i 

(c - 7>Mik/(r # + r„) + £ r ik r kl \ 

*-/ 

+ O' + 1 >M 2 '■/./ + ! “ U| + 7M2 + (c - j)n\]rn 
and thus for / > j 

02) r, v - 10+ 1 >M21 + (c-v')mi Yt r ik r kj\l\k\ + 7M 2 + (1 - r /; - r , t )) 

k-i 

where the sum in equation (12) is defined as zero if j — i - 1. Equation (11) and (12) give an 
easy iterative procedure where the diagonal elements of R are first computed, then the ele¬ 
ments such that j - i - 1 are computed then the elements such that j - i - 2, etc. In this 
manner the exact expression for R is obtained and the solution of equation (9) remains. 

The obvious difficulty with (9) is that for a typical system in which c — 48, the dimension 
of 7" is 2352 x 2352. In order to obtain a solution to (9) it will be reduced to solving a 49 x 49 
system. 


Equation (9) can be rewritten as 


(13) 

Mo + Pi M\ - 0 


(14) 

P*-|A + p k A k + P* + i M k *\~ 0 for k - 1, .. 

. . c - 2 

(15) 

p c - 2 A + Pc-^c-i 4- RM C ) - 6. 



Each of the submatrices in the above equation are of dimension (c + 1) x (c + 1). New 
matrices B k are defined with dimension (2c + 2) x (2c + 2) by 


Bk 


-A k A-' / 

~Mk+\ A' 1 0 f ° r * 


l.c - 2. 
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It now follows from (14) that 

(16) P*> “(**. P*+l>®4 
and thus 

(17) (p 0> p,) - (p,._ 2 , p,._|)fl t _ 2 fl r _ 3 ■■■ B\- 
Combining equations (13) and (15) with (17) results in 

f ^ol 

(18) p,_,(M,_, + RMcW '.DB^i ... fl, w - 0. 

Equation (18) gives a solution for p ( _ t that is unique up to a multiplicative constant. Any such 

solution is obtained and equation (16) is used to obtain p,._ 2 , P,-j.Po The norming 

equation of (9) is then used and the steady state probabilities are determined. 

A unique combination of two previously determined solutions techniques thus yields a 
computational procedure that can be used to solve a class of problems in the telecommunica¬ 
tions field. Typical measures of effectiveness can easily be determined in the same manner as 
in Neuts 16). Thus, the utilization of analytical models for such telecommunication systems is 
now feasible. 
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A NOTE ON THE INFLUENCE OF MISSING OPERATIONS 
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ABSTRACT 

This paper attempts to resolve the existing confusion concerning missing 
operations. Scheduling problems are classified in two groups: (i) null- 
continuous (NCI—comprising the problems where an optimal schedule remains 
optimal on replacement of arbitrarily small processing times (existing opera¬ 
tions) with zeros (missing operations); (ii) null-discontinuous (NDC) — 
comprising those problems which are not null-continuous. 


A "zero processing time" of an operation refers to either of the two following contingen¬ 
cies: (i) An actual operation whose processing time tends to zero. Thus, if r„ denotes the pro¬ 
cessing time of operation 0 7 in job then for a sufficiently small positive number e, scheduling 
problems with r ( , = e and = 0 have the same optimal schedules (algorithms), and e may rea¬ 
sonably be replaced with zero to facilitate calculations, and (ii) A nonexisting (missing) opera¬ 
tion. 


Since an infinitesimal (arbitrarily small) processing time operation (i) has a starting time 
while a missing operation (ii) has not, the two types have a different effect on a scheduling 
problem and must be differentiated to prevent ambiguity. Accordingly, we propose to designate 
an infinitesimal processing time operation as «, and a missing operation as a zero. 

Scheduling problems involving operations with arbitrarily small processing times (existing 
operations where, for all practical purposes, the length may be considered as zero, but have 
starting times) are basically the same as those with the usual strictly-positive processing times 
and do not merit separate consideration. 

Whether a discipline (flowshop or openshop) allows missing operations or not depends on 
its definition, without clear preference for one definition over another. However, for the sake 
of understanding, the definition (whatever it is) must be known and accepted. We propose to 
define flowshop and openshop disciplines allowing missing operations. (Note, [2] and 13] 
define flowshop allowing missing operations while [4], [9] and [14] implicitly assume that 
flowshop does not allow missing operations.) Accordingly, f or 0 in the notation n/m/y/ 5*, 
y € [F,0] comes under that definition and the processing time of the 0, operation of job 
J„ (i — 1,2, .... zr, j — 1,2. m ), is a nonnegalive integer, f (> ^ 0—positive if the opera¬ 

tion exists and zero if it does not. In the case of a flowshop, or an openshop, where all zt jobs 
have m operations—i.e., where the processing times of all mn operations are strictly positive— 

‘The notation is that proposed by Lenstra (9| and Rinnooy Kan 1)4). 
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the constraint r, ; > 0 V ( ij ) is added, and the designation becomes nlmly, t„ > 0/8, 
y ( {f,0|. 

Panwalker, Smith and Woollam {111, demonstrated that the commonly held view that 
there exists an optimal permutation schedule for n/3/F/C mi% and n/m/F\ no-wait/C mas , 
m > 3, turns out to be correct only for positive processing times.* (The latter fact was recog¬ 
nized by Baker (21.) 

1. NULL-CONTINUOUS AND DISCONTINUOUS SCHEDULING PROBLEMS 

The general rule is that optimal schedule remain unchanged under small changes in the 
close neighborhood of the given processing times, providing there is no replacement of e with 
zero (deletion of an infinitesimal operation). Where such replacement is resorted to. some 
problems remain unaffected while in others the optimal schedule undergoes total transforma¬ 
tions; in the latter category, an optimal algorithm and schedule for > 0 V( i,j ) are not 
optimal for t u > 0. Moreover, there are problems for which there exists an efficient optimal 
algorithm (belonging to P) for the first case ( t tj > 0 V (/,/')) while for the second (f„ > 0), the 
problem is NP-Complete. 

DEFINITIONS 

Null-continuous (NC) scheduling problem — one where an optimal schedule remains 
optimal on replacement of arbitrarily small processing times with zeros and vice versa. Thus, 
an optimal algorithm and schedule for /„ > 0 V (ij) are also optimal for t„ ^ 0. 

Null-discontinuous (NDC) scheduling problem — one which is not null-continuous. 
Thus, an NDC problem defines two different problems, one with strictly positive processing 
times (tjj > 0 V (i,j) — missing operations are not allowed) and the other where zero processing 
times are allowed (t 0 ^ 0—missing operations are allowed). 

Informally, a problem where all operations with arbitrarily small processing times for 
every machine can be shifted to the beginning or the end of the schedule without affecting the 
measure of performance or violating the constraints—is NC. To show that a problem is NDC, 
it suffices to work out an example where replacement of arbitrarily small processing time with a 
zero results in a different optimal schedule. 

The following numerical examples of NC and NDC problems demonstrate the fact that an 
NDC problem defines two different problems—one where t 0 > 0 and the other where r„ ^ 0. 

EXAMPLE 1; 

NC instance — 5/2/F/C max + with processing limes as per Table 1. 


Ml I] was brought to our attention by the referee. At the time of writing the first version of the paper, we were 
unaware of 111) and independently reached the same conclusions, 
tThe n/2/ F/C mtt and n/2/F, no-wait/C mix problems are NC and NDC. respectively. 
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TABLE 1 — Processing Times 


Job 

Operatiofbv 
(machine) \ 

12 3 4 5 

1 

2 

4 « e 3 2 

2 2 2 e 4 


The problem is NC since arbitrarily small processing time operations on the first and second 
machines can be shifted to the beginning and the end of the schedule, respectively, without 
violating any constraint or affecting C max . Thus, optima! schedules constructed by Johnson’s 
algorithm (8] remain optimal after replacement of « with zeros. 

NDC Instance— 5/2/F, no-wait/C ma , with processing times again as per Table 1. For 
t 0 > 0 V ( i,j) Gilmore and Gomory’s algorithm [5], 113] yields an optimal schedule as 
presented in Figure 1. The algorithm itself is not affected by replacing e with zeros, but the 
resulting schedule is no longer optimal (Figure 2). 








M, 

1 

El 

Ji 

J4 








m 2 

■ ■ 

m 

J 5 

El 

■ 

m 


0 2 4 6 8 10 11 

Figure I Optimal schedule for 5/2 IF, no-wait. /„ > 0/C ma> , <e > 0) 
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Figure 2. Optimal schedule for 5/2/ F. no-wait/ C mas , (f > 0) 


NP-Completeness of a scheduling problem is proved by reducing a known NP-Complete 
problem to an instance of the problem in question that in many cases contains zero processing 
times. In the light of the preceding discussion only if the problem in question is NC then 
results obtained for /y > 0 V(/j) are also true for t tJ > 0, and vice versa; if NDC, the problem 


VOL 29, NO 3, SEPTEMBER 1982 


NAVAL RESEARCH LOGISTICS QUARTERLY 










538 


\ III I 11/ \M> I \l >IRI 


may refer either to t„ > 0 V(/,/) or /„ ^ 0 with a different approach and specific optimal algo¬ 
rithm and schedule in each case. Thus, before zero processing times are allowed, the problem 
in question must be characterized as NC or NDC 

2. CLASSIFICATION AND COMPLEXITY OF NC AND NDC PROBLEMS 

A first attempt has been made to classify some of the scheduling problems into NC and 
NDC, and to examine the complexity of each problem (the results are obtained in Table 2). 
As explained previously, each NDC problem defines two different problems, one with strictly 
positive processing limes (r„ > 0 V(/',/)) and the other with nonnegative ones, (f„ > 0). The 
fact that for a particular problem, the case where /„ > 0 V (/,_/) is a special case of t,, ^ 0 and. 
thus, its complexity is lower than or equal to that of the latter (theoretically speaking, all com¬ 
binations are possible) is demonstrated in Table 2. 


TABLE 2 — Classification and Complexity of NC and NDC Problems 



Problem 

Complexity 


r„ > 0 

> o 

u 

n 

z 

n/2/F, no-wait/C max * 

n/3/A, no-wait/C max 

n/4/F, no-wait/C max 

n/2/F, r, > 0/C max 

n/2/F. tree/ C max t 

zr/2//vtree7C max § 

n/3/F/C max 

n/ 2/0, no-wait/C max * 

n/2/J, n, = 2, no-wait/C max * 

n/2/FfLC, 

n/2/F, no-wait/IC, 

n/ 2/0, no-wait/IC, 

n/2/0/IC, 

n/2/J/ni = 2, no-wait/IC, 
n/l/seq.dep./ C max 
(Travelling salesman) 

0</; 2 >, [5] 

Unary NP-Complete, [12] 
Binary NP-Complete, [10] 
Binary NP-Complete, [10] 
0(n log n), [16] 

Unary NP-Complete, [10] 
Unary NP-Complete, [15] 
Unary NP-Complete. [15] 
Unary NP-Complete, [4] 

Of 

0 

0 

0 

Unary NP-Complete, [14] 

Unary NP-Complete, [151 
Unary NP-Complete, [4] 
Unary NP-Complete, [121 
Binary NP-Complete, [10] 
Binary NP-Complete. [10] 

‘> 

Unary NP-Complete. [4] 
Unary NP-Complete, [15] 
Unary NP-Complete. [15] 
Unary NP-Complete, [4] 
Unary NP-Complete, [4] 
Unary NP-Complete, [1] 
Unary NP-Complete, [1] 
Unary NP-Complete. [4] 
Unary NP-Complete.[14] 

u 

z 

n/2/ F/C max 
n/2/J, n, < 2/ C max 
n/m/ 0/8 

0(n log n), [8] 

0(n log n), [7] 

Depends on the number of machines (m) and the regular 
measure of performance (8). 


'Elowshop. jobshop and openshop disciplines are defined allowing zero processing limes 
+ Prize carrying open problem. [6|, (9|. |I0| 

*lree—Tree precedence relations where job 7, may starl only after job 7, has been completed 

§trce’ - Tree p' cedence relations where, on each machine, job 7, may start only after 7 on this machine 

has been completed 
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